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Abstract 

Granular fluids consist of collections of activated mesoscopic or macroscopic particles (e.g., pow- 
ders or grains) whose flows often appear similar to those of normal fluids. To explore the qualitative 
and quantitative description of these flows an idealized model for such fluids, a system of smooth 
inelastic hard spheres, is considered. The single feature distinguishing granular and normal fluids 
being explored in this way is the inelasticity of collisions. The dominant differences observed in 
real granular fluids are indeed captured by this feature. Following a brief introductory descrip- 
tion of real granular fluids and motivation for the idealized model, the elements of nonequilibrium 
statistical mechanics are recalled (observables, states, and their dynamics). Peculiarities of the 
hard sphere interactions are developed in detail. The exact microscopic balance equations for the 
number, energy, and momentum densities are derived and their averages described as the origin 
for a possible macroscopic continuum mechanics description. This formally exact analysis leads 
to closed, macroscopic hydrodynamic equations through the notion of a "normal" state. This 
concept is introduced and the Navier-Stokes constitutive equations are derived, with associated 
Green-Kubo expressions for the transport coefficients. A parallel description of granular gases is 
described in the context of kinetic theory, and the Boltzmann limit is identified critically. The 
construction of the "normal" solution to the kinetic equation is outlined, and Navier-Stokes order 
hydrodynamic equations are re-derived for a low density granular gas. 

These are notes prepared as the basis for six lectures at the Second Warsaw School on Statistical 
Physics held in Kazimierz, Poland, June 2007. 

PACS numbers: 



1 



I. INTRODUCTION 



Granular fluids are ubiquitous in nature l|, l2|. To illustrate with the familiar first, go to 
your kitchen and take out the mustard seeds, pepper corns, salt, and rice. Put 100 grains 
of each in small jars and tumble these jars in various directions. The gravitational field 
induces temporary flows of groups of grains in each case. While the flows appear to be 
locally convective it is clear that each grain's motion is more complex due to interactions 
with other grains. This is qualitatively similar to macroscopic flows of normal fluids (those 
composed of atoms or molecules) in which the coarse grained convection is the collective 
effect of complex atomic collisional motion. The objective here is to explore to what extent 



the well developed methods to describe the macroscopic dynamics of normal fluids [3| can 
be extended to the flow of such compositions of grains 1^. Clearly, there are signiflcant 
differences in the physical systems. It will be seen that many properties of interest are 
insensitive to some of the most obvious differences, and that others can be incorporated in 
realistic models for the system of grains. 

Normal fluids incorporate a wide range of different systems. Simple atomic fluids are well 
represented by spherically symmetric central point forces between the atomic constituents, 
and Navier-Stokes hydrodynamics (deflned precisely below) applies to most macroscopic 
nonequilibrium states of interest in this case. The interactions in low molecular weight 
fluids no longer have spherical symmetry, but at the macroscopic level this is a quantitative 
rather than a qualitative effect (e.g., only the values of transport coefficients in the Navier- 
Stokes description change). Mixtures of these types of fluids also have the same qualitative 
macroscopic behavior. On the other hand complex fluids, such as those composed of high 
molecular weight (polymers), generally exhibit quite different macroscopic behavior which 
can depend on both the system and the class of macroscopic states considered. 

Some of this diversity is evident for systems of grains as well. The mustard seeds are 
more mono disperse and smooth than the pepper corns, but both are roughly spherical. 
On the other hand rice is asymmetric to differing extents (e.g., large aspect ratio in China, 
small in Spain). As with normal fluids, their macroscopic flows are quite similar. Salt has 
irregular shape and furthermore is complicated by being hydrophilic - any moisture in the 
air changes dramatically the interaction between grains. This can change its motion as well, 
as everyone knows about salt shakers at the seaside. Thus, the medium in which the grains 



2 



move can be important. Structurally complex grains such as collections of filaments will 
not be considered here. Still, structurally simple granular systems can behave as complex 
fluids in many nonequilibrium states of interest. This is an important difference between 
normal fluids and granular systems that provides some of the most difficult challenges and 
also some of the most interesting opportunities. 

Beyond the kitchen, granular systems include objects of interest to the phamacutical 
industry (pills and associated powders generated in their production), the agricultural in- 
dustry (storage and transport of edible grains), geology (rock and snow avalanches), and 
extra-terrestrial systems (Saturn's rings, regolith on Mars). Two general classes of states 
for these granular systems are distinguished, compact and activated. Unshaken, the rice in 
the jar appears at rest. In fact, each grain has some kinetic energy due to the temperature 
of the room. However, on Earth the gravitational potential energy relative to the bottom of 
the container for heights h greater than the dimension of the grain is much larger than this 
energy of motion, due to the large mass of typical grains. Thus, they pack at the bottom 
of the container. Questions about the possible distribution of packing configurations con- 
stitutes a field of current activity [s], and is of central practical interest for the storage of 
grains. For example, granular storage in silos commonly leads to explosions whose mitiga- 
tion requires knowledge of the distribution of forces in the packed grains. Here, attention is 
limited to the second class of activated grains. Work is done on the system (e.g., shaking) 
to provide a kinetic energy greater than that required to overcome the compactification by 
gravity. In addition, when the number of grains is large the activation typically induces 
an apparent random component to the granular motion due to frequent collisions among 
the grains. Systems of activated, collisional grains constitute the qualitative definition of 
granular fluids for the discussion here. 

Phenomenologically, Navier-Stokes hydrodynamics has been applied to a wide class of 
granular flows in practice with qualitative success in many cases. To what extent can such a 
description be justified from a more fundamental basis? If justified, how can the parameters 
of this description (equation of state, energy loss function, transport coefficients) be given 
fundamental definitions as functions of the state conditions? If justified, what is the context 
and limitations of this description? These are the issues that are addressed here. As with 
the history of progress on these same questions for normal fluids, initial attention is focused 
on an idealized fluid and states where conceptual problems can be isolated and addressed 
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cleanly. Steps to describe more realistic granular fluids can then proceed with increased 
confidence and guidance. In the next section this idealized granular fluid is deflned and the 
overview of its exploration is summarized. 

II. OVERVIEW OF THE PRESENTATION 

The mustard seeds are hard, spherical, and monodisperse. This system can be refined for 
experiments using carefully machined spheres of glass or metal with empirically determined 
binary collision properties. Although hard, the binary collisions are quite complex in detail 
since each is composed of a large number of composite molecules. While in contact, the 
shape of each grain is distorted and energy is redistributed between the kinetic energy of 
their centers of mass and that of their internal degrees of freedom. On separation, they again 
move freely but with some energy lost to internal and rotational degrees of freedom. The 
"hard" interaction means that the contact time is short, so effectively the motion is that of 
free streaming punctuated by velocity changes on collision with a consequent loss of energy 
in each case. This suggests the idealization of a system of hard spheres (zero collision time) 
with inelastic collisions. A further simplification is the neglect of transfer between kinetic 
and rotational motion. This idealized fiuid then consists of a system of smooth, inelastic 
hard spheres. It is similar to the idealized hard sphere fiuid for normal atomic systems, with 
the sole difference being an energy loss on pair collisions. 

In the next section the context of this ideal granular fluid is elaborated further to clarify 
the extent to which its properties should represent both qualitatively and quantitatively 
many properties of real granular fluids. The elements of nonequilibrium statistical mechanics 
are recalled briefly and applied to this system of hard, inelastic, smooth spheres 
. Since the forces are singular, the usual form for the generator of dynamics for piecewise 
continuous forces no longer applies and the necessary changes are described in Appendix A. 
As a consequence of these changes the generators in different representations for the time 
dependence are different (e.g., that for the Liouville equation and that for the observables) . 
Three different generators are identifled in Appendix A. The average energy for an isolated 
system is shown to decrease monotonically due to the loss of energy on each binary collision. 
Consequently, there is no equilibrium Gibbs solution to the Liouville equation. Instead, 
it supports a "homogeneous cooling state" (HCS) in which the dynamics of the energy 
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loss appears only through a scaling of the velocities (hence the terminology "cooling") and 
associated normalization factors. Empirically (i.e., in molecular dynamics simulations) it 
is found that a wide class of homogeneous states approach the HCS after a few collisions 
per particle Consequently, the rapid velocity relaxation in normal fluids leading to the 
Gibbs distribution has its counterpart in this granular fluid in the approach on a similar time 
scale to the HCS distribution. The scaling form of the HCS suggests a related dimensionless 
representation for the Liouville equation in a more general context. 

As described above, a main target of this investigation is the possibility of a macro- 
scopic hydrodynamics. This refers to closed equations for the hydrodynamic observables: 
number density, energy density, and momentum density. As a fundamental starting point, 
the exact microscopic phase functions corresponding to these fields are identified and the 
balance equations for them are obtained in Appendix B. For a normal fluid, these are the 
local microscopic conservation laws. Section HVl discusses the average of these equations, the 
macroscopic balance equations, which is the framework in which the corresponding hydrody- 
namic equations can exist. Contributions to the fluxes in these equations due to convection 
are identified from a local Galilean transformation p^. The remaining contributions are 
due to "dissipative" exchanges of density, energy, and momentum for fluid cell at rest. The 
derivation of explicit expressions for this remainder, called "constitutive equations", is the 
central problem for obtaining hydrodynamic equations. A change of variables from density, 
energy density, and momentum density to density, temperature, and flow velocity is defined 
and introduced. Finally, the macroscopic state corresponding to the HCS is identified from 
these equations. 

A more precise definition of hydrodynamics is considered in Section |Vl where the notion 
of "normal" states is introduced and motivated. It is shown how this concept leads to 
constitutive equations which a functionals of the hydrodynamic fields, giving in principle 
a set of five closed equations for the determination of these fields. The special case of 
weakly inhomogeneous states is considered in Section I VI I as an example, leading to the 
phenomenological Navier-Stokes equations for a granular fluid. Two important differences 
between normal and granular fluids are noted, long wavelength instability and the existence 
of new stationary states. 

The formal construction of a normal solution to the Liouville equation is considered in 
Section rVIII for weakly inhomogeneous states by an expansion in the local spatial gradients. 
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The leading order reference state is a local HCS distribution, similar to the local equilib- 
rium Gibbs distribution for normal fluids. The constitutive equations at this approximation 
provide definitions for the hydrostatic pressure and the local cooling rate in each fiuid cell. 
The formal solution to the Liouville equation including first order in the gradients provides 
constitutive equations characterized by the transport coefficients for a granular fiuids. This 

Tom a solution to the Liouville equation for small initial 
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12| . Expressions for these transport coefficients, the 



extends recent results obtained 
perturbations (linear response) 
Green-Kubo formulas for granular fiuids, also are discussed in Section IVIII The specific 
case of shear viscosity is explored further for illustration. The simplest test of these formal 
developments is provided by the dynamics of a single impurity in a granular fiuid in its HCS, 
corresponding to impurity motion in an equilibrium fiuid. The expected hydrodynamics in 
this case is simple diffusion of the particle's probability density. The Green-Kubo expres- 
sion for the diffusion coefficient is noted, and a practical short time approximation for the 
associated velocity autocorrelation function is demonstrated [l3^. 

A different approach to the determination of constitutive equations is provided by ki- 
netic theory, a representation in terms of the reduced single particle distribution function. 
This approach is described in Section IVlllI where the objects of interest are shown to be 
determined exactly by the single particle distribution function, and the hierarchy of equa- 
tions determining all reduced distribution functions is derived from the Liouville equation. 
A formal solution to the hierarchy is obtained systematically at low density, resulting in 
a Boltzmann description for a granular gas. The derivation of hydrodynamics from this 
kinetic theory is outlined and compared to the general fluid results. 

Finally, some additional perspective on these results is given in the Discussion Section. 
The analysis here suggests that hydrodynamics is a useful concept for granular fluids, but its 
context and limitations have not been explored in any detail. Appropriate space and time 
scales for hydrodynamics in general and limitations of the Navier-Stokes approximation in 
particular are discussed. 

As these notes are prepared for a series of lectures the material is focused on the research of 
the author and his colleagues. The references quoted are heavily weighted in that direction 
as well. Apologies are offered at the outset for the many excellent contributions to this 
subject matter not given explicit recognition. 
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III. IDEALIZED GRANULAR FLUID AND STATISTICAL MECHANICS 
A. Hard sphere idealization 

Consider a fluid comprised of mono- disperse, spherical particles with pairwise additive 
central interactions (i.e., smooth particles with no tangential momentum transfer). Also, 
assume that Newton's third law holds so that momentum is conserved between colliding 
pairs. More general cases of mixtures, non-central, or many-body forces can be incorporated 
with greater complexity but no signiflcant conceptual changes. The case of dissipative " soft" 
spheres is considered flrst. The pair force is assumed to be piecewise continuous for relative 
distance between particles r and s, Qrs = |qr ~ Qsl < and vanishing for g^s > cr- As a 
central force it is directed along the line of centers and therefore has the form 



Here grs = — is the velocity of approach or separation. From spherical symmetry the 
magnitude of the force can depend only on the scalars g^s? Ors, firs' grs- The functional form of 
this force is such as to describe the two physical effects of repulsion and dissipation during 
the deformation. For example, it could be the superposition of an conservative repulsive 
elastic force plus a drag force proportional to the normal component of the relative velocity. 
The amount of deformation is then d ~ \/eJk where e is the energy per particle and k is the 
elastic constant. The conditions of interest here are such that d/cr << 1 and Tc/rm << 1, 
where ~ d^fmje is the average contact time and ~ rT^I'^ ^fmje is the mean free time 

1 /3 

between collisions. Consequently, {Tc/Tm) ~ (ncr^) {d/(y) < (d/a) since na^ < 1 for fluid 
states. These are rough estimates, but they show that the controllable parameters e and k 
admit conditions where the particles behave as hard objects. A force law / (g^s) ~ i'^/^)"^ 
describes well the repulsive forces for simple atomic fluids, with n > 12, and its properties 
are known to be accurately represented by those for hard spheres n — » oo [3]. A similar 
limit for the conservative part of the deformation potential can be expected as well. 

It is necessary also to retain the fact that there is a total energy loss during the contact 
time in this hard sphere limit. This is done by requiring that the collisions are inelastic. 
In detail, the total momentum for the colliding pair is conserved and the relative velocity 
changes according to 




(3.1) 



g. 



'rs 



<rs 



Qrs ' grs ) ^ ) 



(3.2) 
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with the corresponding energy change 



A Qm (f^ + = Im {gll - gl^) = -Kn{l-a^ (q,., ■ g^,)) (q^, ■ g,.,)^ . (3.3) 

The scalar parameter a (q^s ■ grs) is the restitution coefficient and takes on the values < 
a < 1, with a = 1 corresponding to elastic collisions. A detailed correspondence between 
inelastic hard spheres and soft dissipative spheres requires that the restitution coefficient 
should depend on the normal component of the relative velocity iJ]. In particular, it is 
found that a (q^s ■ grs) ^ 1 as q^s ■ grs 0. However, if the conditions studied avoid this 
limit (e.g., continual energy input), then it is possible to consider the further idealization 
of an average energy loss per collision characterized by a constant restitution coefficient 
a {f\rs ■ grs) — * «• This will be the case studied in all of the following. 

There is a mathematical price to be paid for this idealization of smooth, inelastic, hard 
spheres. The forces are singular at contact, and therefore the usual description of the 
dynamics from Hamilton's equations must be modified fiol . [l5 . [igI . \\\ . The modification 
corresponds to replacing the effect of the continuous force by a binary scattering operator 
that generates the instantaneous momentum change of the pair on contact. The form of 
this operator is derived in Appendix A. This complication is common to the representation 
by hard spheres of both normal and granular fluids. In the latter case there is another, 
perhaps unexpected, effect of the hard collisions plus dissipation called inelastic collapse 



18( . To understand this, consider a perfectly elastic hard ball dropped from a height h on 



a horizontal table in the gravitational field. Without interference it bounces indefinitely and 
it is never in contact with the table over any finite time interval. In contrast, the inelastic 
hard sphere will undergo an infinite number of collisions in a finite time r and come to rest 
in contact with the table. A similar effect occurs in the inelastic hard sphere granular fluid 
where groups of particles can cluster with increasing numbers of collisions among them in a 
finite time interval. The issue of inelastic collapse will be revisited at the end of this section. 



B. Overview of nonequilibrium statistical mechanics 

For a given physical system its description via statistical mechanics consists of states 
(represented by distribution functions over its phase space), observables (phase space func- 
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tions), and expectation values (averages of the observables over the states). Its evolution 
in time is given by a mapping of either the states or the observables in phase space. The 
system of interest here is a one component fluid of identical smooth, inelastic, hard 
spheres (mass m, diameter a). The state of the system at time t = is completely char- 
acterized by the positions and velocities of all particles, represented by a point in a 6A^ 
dimensional phase space Fq = {qi(0), . . . , qAr(O), vi(0), . . . , Viv(O)}. The dynamics consists 
of straight line motion along the direction of the velocity at time t (free streaming), until 
any pair of particles, say is in contact. The relative velocity gjj = Vj — Vj of that pair 
changes instantaneously according to a given collision rule (13. 2p while their total momentum 
is unchanged. Subsequently, all particles continue to stream freely until another pair is at 
contact, and the coUisional change is repeated for that pair. In this way a unique trajectory 
Tt = {qi(t), . . . , qAr(t), vi(t), . . . , VAr(t)} is generated in the phase space for t > 0, where 
the configurational degrees of freedom change continuously while those for the velocities are 
piecewise constants. 

The statistical mechanics for a fluid of inelastic hard spheres has been described elsewhere 



0, 0, y. 
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20| . It is comprised of the dynamics just described, a macrostate specified 



in terms of a probability density p{T), and a set of observables denoted by A{r). The 
expectation value for an observable at time t > for a state p(r) given at t = is defined 
by 

(A(t);0)= I dTp{T)A{T,) (3.4) 

where A(t) = A(Tt), and Tt = {qi{t), . . . , qAr(t), vi(t), . . . , VAr(t)} is the phase point evolved 
to time t from F = Tt=Q. The dynamics can be represented in terms of a generator L defined 
by 

{A{t);0) = j dV p{T)e'^A{T). (3.5) 

For continuous potentials the generator is easily recognized from Hamilton's equation as a 
Poisson bracket operation with the corresponding Hamiltonian. However, its identification 
for the discontinuous hard sphere potential is less direct [lO, 15, 16, 17|. There are two 
components to the generator, corresponding to the two steps of free streaming and velocity 
changes at contact. The first part is the same as for continuous potentials while the second 
part replaces the contribution from the singular force by a "binary collision operator" T{i,j) 
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for each pair J 21 1 



N , N N 



i^ = 5^v,-V, + ^5^X^mj). (3.6) 

1=1 1=1 jj^i 

The binary coUision operator T{i,j) for normal fluids is identified directly from the Poisson 
bracket of Hamilton's equations 

T(z, j) ^ % = m-'F{q,,) ■ (V,^ - V.J . (3.7) 

where F{qij) is a conservative force. For hard spheres, the position variables are still contin- 
uous functions of time but the momenta are piecewise constant (in the absence of external 
forces) and discontinuous. The form for T{i,j) in this case is obtained in Appendix A with 
the result 

TiiJ) = Qi-gij ■ ■ q.ij\5iqij - (T){bij - 1), (3.8) 

where 9 is the Heaviside step function, and bij is a substitution operator which changes the 
relative velocity gij into its scattered value g^^, given by Eq. (13.21) 

%A(g,,)=A(g^)- (3.9) 

The theta function and delta function in (13.81) assure that a collision takes place, i.e. the 
pair is at contact and directed toward each other. 

An alternative equivalent representation of the dynamics is obtained by transfering the 
dynamics from the observable A{r) to the state p(r) by the definition 

J rfrp(r)e*^A(r) = ^ (e"*^p(r)) a{t). (3.10) 

The representation in terms of a dynamical state is referred to as Liouville dynamics. The 
probability density p(r) must vanish for all configurations of overlapping hard spheres, so 
the domain of integration on the left side of fl3.10p is effectively restricted to non-overlapping 
configurations. Thus the generator L is used always in that context. However, the right 
side of (13.101) no longer has that restriction and consequently the generator for Liouville 
dynamics is not the same as that for observables (as in the case of continuous potentials). 
Instead, direct analysis of (I3.10p leads to the result (see Appendix |A]) 

N ^ N N 

I = ^ V, . V. - -J2Y.Tihj), (3.11) 

1=1 1=1 j^i 
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with the new binary colhsion operator 

T_{i,j) = 6{qij - a)\gij ■ q^j\{&{gij ■ q,ij)a^\^ - Q{-gij ■ %)). (3.12) 

Here b^j^ is the inverse of the operator bij in fl3.9p . 

Next consider time correlation functions for two observables A and B 



{A{t)B;0) 



J dT (e*^A(r)) p{T)B{T) = J drA(r)e-*^ {p{T)B{T)) . (3.13) 
A third generator is defined for reversed dynamics along the trajectory by 

e-*^ (p(r)S(r)) = (e-*^p(r)) (e-*^-i?(r)) . (3.14) 

This gives 

{A{t)B-0) = J dTA{T)e-'^ {piT)B(T)) = j dV p{T ,t)A{T)e-'^- B{T) 

= {AB{-t);t). (3.15) 

where the reversed dynamics of B{—t) has been defined by 

B{-t) = e-'^-B{T). (3.16) 

The new generator is identified in Appendix |X] as 

N ^ N N 

L_ = 5^ V. ■ V,, - - 5^ 5^ T_(z, j) (3.17) 

with 



2 

4=1 i=l jjii 



T_{i,i) = 6{qij - cr)Q{gij ■ qij)\gij ■ (6-/ - l) (3.18) 

In summary, the problems presented by the singular forces for a fluid of hard spheres are 
resolved if Hamilton's equations for observables are replaced by 

{dt - L) A{T, t) = 0, {dt + L_) A{T, -t) = (3.19) 

and the Liouville equation for probability densities is replaced by 

{dt + L}piT,t) = 0, (3.20) 

for t > 0, with the respective generators given by (13. 6p . (13.171) . and (13. lip . The three gener- 
ators are all different, so some care must be used to apply them under the correct conditions 
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of their definitions. Note that the forms of the generators L and L_, and corresponding 
binary coUision operators T{i,j) and T^{i,j), do not depend on the details of the colhsion 
rule defining the operator bij; the result applies for both elastic and inelastic collisions 21 1. 
In contrast, the generator for Liouville dynamics is obtained by a change of variables that 
introduces the Jacobian of the transformation between the variables gij and bijgij. Hence it 
depends explicitly on the collision rule and the restitution coefficient a. 

For normal fluids the probability of a configuration with any pair of particles in contact 
and at rest has vanishing measure. Then the above description of trajectories as sequences 
of pair collisions among the particles is adequate. As noted above, in granular fiuids the 
phenomenon of inelastic collapse admits the possibility of evolution to a state where clusters 
of particles are in contact and at rest. It would appear that the collision of another particle 
with that cluster be described would require a more complex generator for the dynamics 
than that considered here. Consider the simplest case of a pair in contact and at rest, with 
a third particle incident on one of the two. The generator here would transfer momentum 
to one of the pair, as though it were isolated. Then, it would no longer be at rest with 
respect to the other member of the pair and a second instantaneous momentum transfer 
would occur to the second member. As a result, both particles originally at contact and in 
relative rest would experience relative motion and all three particles would separate. In fact, 
this is the correct dynamics for real particles and it avoids the difficulty of indeterminate 
dynamics for collisions with clusters. Mathematically, therefore, there appears to be no 
difficulty for the dynamics generated here due to inelastic collapses. In practice, however, 
this can be difficult to simulate as the effective time between such collisions can be very 
small and require following a large number of collisions. 



C. Homogeneous cooling state 

In the absence of external forces there is special solution to the Liouville equation for 
normal fluids: the stationary, homogenous (translationally invariant) equilibrium solution, 
Pg. For the isolated system considered here this is a probability density with sharply defined 
total energy, total momentum, and number of particles. It is "universal" in the sense that 
most other homogeneous initial preparations rapidly approach this stationary equilibrium 
solution, on a time scale of the order of several collisions per particle. This is the collisional 
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velocity relaxation to Maxwellian distributions for each velocity degree of freedom. Granular 
fluids are different in the sense that the Liouville equation for an isolated system has no 
stationary solution. This is due to the loss of energy on each inelastic collision, such that 
the total energy decreases monotonically E{t) < E{0). Nevertheless, it appears there is a 
universal homogeneous solution whose time dependence occurs entirely through a scaling 



of the velocities for each particle 
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Here q^^ = — q^, and \Jh is an overall constant velocity of the system. This velocity can 
be removed by a Galilean transformation but it is useful to retain it for the generalization 
to a corresponding local form defined below. Also, / is an arbitrary constant characteristic 
length and Vh{t) is a "thermal" velocity defined in terms of the energy per particle 



2,,, 2n{t) ^ _ 2Eh{t) 2 



Kit) = Ut) ^ ^ = ^ ldT\y^ -mvt 1 p, (r; t) (3.22) 

The second equality defines an associated "temperature" for the system, so that Vhit) can 
be interpreted as the average thermal speed. This temperature definition agrees with that 
introduced below for one of the hydrodynamic fields. Then, (F; t) = (F; T{t)) is an 
example of a "normal" state, also defined below, whose time dependence occurs entirely 
through the hydrodynamic fields. 

The specific form for p^ is determined by the requirement that Pf^ should be a solution 
to the Liouville equation f l3.20p 

{dtT,)dT,p, + Lp, = 0, (3.23) 

or 

£t,P, = 0, Zt = -C,(T)T9t + X (3.24) 
The " cooling rate" C/i has been introduced by the definition 

CAThit))^-T,\mUt). (3.25) 

The scaling property (13.211) gives an equivalent alternative form 

Zp, = 0, (3.26) 
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where the operator £ is defined by 

- r ^ 

CX = ^J2^-r- [(v. - U,) X] + LX. (3.27) 

r=l 

An exphcit expression for the coohng rate in terms of follows from differentiation of fl3.22p 
with respect to time and using the Liouville equation to get 

= -iN-l)-^JdTp,{T;t)T{z,j){v^ + v^). (3.28) 

Recalling the explicit form for T{i,j) in (13.81) . the cooling rate is seen to be proportional to 
the energy changes on binary collisions given by (13.31) . The cooling rate then simplifies to 

Nm f 

C. = (1 - «')Y^jr^ y rfrp,(r,t)(gi2-qi2)'e(gi2 ■ qi2)5(gi2 - o). (3.29) 

Note that since (r;t) depends on time only through Thit) the cooling rate also has this 
property, as the notation in (13.251) implies. 

In summary, the special HCS solution to the Liouville equation p^^ (F; t) is defined by 
( ]3.26p together with ( ]3.29p which is a linear functional of pi^{T\t). In this representation, 
the linear Liouville equation in terms of r,t becomes a nonlinear equation parameterized 
by Th{t). In principle, the solution is determined in two steps. First, (13.261) and (I3.29P 
are solved as a function of Thit). Second Th(t) is determined from (13.251) and substituted 
into the solution found in the first step. The time dependence from this second step can 
be determined quite generally from the scaling property of his dependence on Th{t) can be 
made explicit through the scaling form of (I3.2ip which leads to 

c. (nit)) = (3.30) 

where is a constant 

a = (1 - «')f / drpiir)ig*u-^i2mgu ■ ^12)^2 - j)- (3.31) 

The dimensionless variables for the integration here are given by 

F*^{qt,..,q^,vt..,v^}, < = ^, v,* = (3.32) 

I Vh[t) 
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Thus fl3.30p exposes the exphcit temperature dependence of (Thit)) and (13.250 becomes 



which can be integrated to give 



n{t) = T,(o) ( 1 + ^-^t 



2m 



-2 



(3.33) 



(3.34) 



The initial time has been chosen at t = without loss of generality. The long time cooling 
is seen to be algebraic and universal (independent of initial conditions). Thus, p^j (X',t) is a 
universal function of Th{t) which itself becomes universal at long times. 



D. Dimensionless representations 

These last considerations suggest that the mathematics may be simpler and the physics 
setter exposed by using a representation in terms of the dimensionless variables (I3.32p 
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201 1 . The associated dimensionless time is defined through the differential form 



ds = "^dt. (3.35) 



This can be integrated using ( 13.340 to give 



sM = ^^4"-^^]- (3-36) 



c;; V 2£ 

The parameter s(t, 0) is a measure of the average number of collisions per particle in the 
interval (0,t). In terms of this parameter the cooling of (13.340 becomes exponential 



The dimensionless form of (13.261) is 



e-^'^^ (3.37) 



Cpl = Q, (3.38) 

with 



£*X = ^ 5^ Vv; ■ [v;X] + L*X. (3.39) 



N 



2 



r=l 



This equation is supplemented by the definition of C,\ in terms of p*i^ in (13.310 . There is no 
longer any time dependence, no dependence on Th{t), only a parameterization of the solution 
by the scalar (l^. Further elaboration on this is given below. 

15 



Now return to the more general solutions to the Liouville equation given by (13.200 . Each 
solution will be characterized by an initial total energy, momentum, and particle number. 
For these global parameters there is an associated and Th{t) defined. The equation for 
general solutions can be expressed in terms of the same dimensionless variables fl3.32p and 



(I3.35P for the associated HCS. The result is [20 1 



ds + £* j p* = 0, p {T- 1) = {IvH (t))-^V (r*; s) . (3.40) 

The dimensionless generator for the dynamics, £*, is the same as that given in fl3.39p . This 
result is similar in form to the original representation of the Liouville equation, except that 
now the time is measured in terms of the average number of collisions and the generator for 
the dynamics has an additional contribution compensating for collisional cooling as it would 
occur in the corresponding HCS. Note that in this representation p*^ is a stationary solution 
to the Liouville equation fl3.40p . The differences between normal and granular fluids have 
been somewhat mitigated by this dimensionless representation. For example, notions of 
"approach to equilibrium" can be translated into "approach to the HCS", and the universal 
features of the equilibrium state can be translated to those of the HCS. MD simulations 
suggest that these comparisons are useful in the sense that very different homogeneous 
initial conditions approach the HCS on a time scale of several collisions per particle [Sj] . The 
corresponding dimensionless representations for observables corresponding to (13. 190 are 

(9, -/:*)A*(r*;s) =0, {ds + CL) A* {V*--s) = (3.41) 
£*X = J]v; ■ Vv;X + L*X, C*_X = ^^wl-V^.X + L*_X (3.42) 

r=l r=l 

IV. MACROSCOPIC BALANCE EQUATIONS 

For a simple one component fluid the relevant macroscopic variables are the average 
number density n, energy density e, and momentum density g 

n{v,t) = {n{v)-t), e(r,t) = (e(r);t), g(r, t) = (g (r) ; t) . (4.1) 

The brackets denotes an average over the ensemble p(r,t) as indicated on the right side of 

dsn 

(X;t)^ j dV (e-*Vr))x(r). (4.2) 
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This set will be referred to as the hydro dynamic fields, as it is their macroscopic dynamics 
which is the candidate for a closed description on some appropriate length and time scale. 
The specific forms of the phase functions are given in Appendix [B] (a caret has been intro- 
duced to distinguish the microscopic and macroscopic fields). Also in that Appendix the 
exact microscopic balance equations for the phase functions are derived. Their averages give 
the corresponding macroscopic balance equations 

atn(r,t) + m-Vr-g(r,t) = (4.3) 

dte{r,t) + Vr ■ {s{r);t) = {w{r)-t) (4.4) 

dtgair,t) + Vr, {hapir); t) = 0. (4.5) 

Here (s(r,i(:);0) is the average energy fiux, (w(r,t);0) is the average energy loss function, 
and {hai3{r,t);0) is the average momentum fiux. Again, the phase functions s(r,t), w(r,t), 
and hai3{r,t) are given exphcitly in Appendix [B] [3]. These exact equations are the starting 
point for investigating the possibility of a hydrodynamic description for a granular fiuid. 

To further simplify the balance equations it is useful to define the average fiow velocity 
according to 

g(r,t) = n{r, t)mU(r, t). (4.6) 

Then as shown in Appendix O the purely convective parts of the energy and fiuxes can be 
extracted by a local Galilean transformation 

e(r, t) = e'(r, t) + ^mn(r, t)U\r, t) (4.7) 

Mv)-t) = {s'^{v)- 1) + ?7,(r, t) (^e'(r, t)+^mn(r, t)U\v, t)^ + {h'^f,{v)- 1) Up{v, t) (4.8) 

( V(r); = (^a/3(r); t) + mn(r, t)f/„(r, t)Up{v, t) (4.9) 

The phase functions with a prime denote the same function evaluated at Vj ^ = Vj — 
U(r,t). Therefore, e'(r, t) is the rest frame "thermal energy", (s'(r);i(:) is the rest frame 
"heat fiux", and {h'^p{v)]t^ is the rest frame "pressure tensor". This terminology does not 
imply any thermodynamic implications, however, and it is useful to retain the names for 
comparison with normal fiuid forms. In this spirit, the following notation is introduced 

e'(r, t) ^ ^n(r, t)T(r, t), C(r, t) ^ -^-^-1^^^ {w{r); t) , (4.10) 
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(s'(r); t) ^ q(r, t), {h'^^{r)- 1) ^ P,^(r,t) (4.11) 

Equation fl4.10p defines the "granular temperature". As a definition, T{r,t) is meaningful 
for any state of the system as a measure of the local kinetic energy e'{r,t) - it simply 
constitutes a change of variables from the pair (n, e',U) to {n,T,\J) as the hydrodynamic 
fields of interest. However, its relationship to any given measuring device ("thermometer") 
must be considered with care. Similarly, (14. lip gives a microscopic definition for the heat 
flux q(r, t) and pressure tensor PQ,^(r,t). The interpretation of C(r,t) as a "cooling rate" 
appears in (14.130 directly below. 

In terms of these new variables the macroscopic balance equations become 

Dtn{v, t) + n(r,t) Vr ■ U(r, t) = (4.12) 

(A + C(r, t)) T(r, t) + (P„/3(r, t)d^Up{r, t) + V • q(r, t)) = 0, (4.13) 

DtU^ir, t) + (mn(r, t))-^dpPo^p{v, t) = 0, (4.14) 

where Dt = dt + \J ■ V is the material derivative. These macroscopic balance equations are 
still exact. They have the same form as for a normal fluid, with the only change being the 
presence of the cooling rate. These are not a closed set of equations for n, T, and U since 
the heat flux, pressure tensor, and cooling rate have not been explicitly represented. Clearly, 
however, at this point these equations apply equally well to both granular and normal fluids 
under the most general fluid state conditions. 



V. "NORMAL" STATES, CONSTITUTIVE RELATIONS, AND HYDRODY- 
NAMICS 

The most general notion of a hydrodynamic description is a closed set of equations for the 
hydrodynamic fields, -v^ {n,T,\J} . This follows from the exact macroscopic balance 
equations if the cooling rate, heat flux, and pressure tensor can be represented as functional 
of these fields 

C(r, t) ^ C(r, t I {y^}), q(r, t) q(r, t \ {yj), P^pir, t) P,^(r, t \ {y^}) (5.1) 

These are known as constitutive relations. A comment on notation is appropriate at this 
point: / (r,t, {ya (r,^)}) denotes a function of T,t and the fields {ya (r, i)} at the point r, 



18 



while /(r, t | {Va}) denotes a function of T,t and a functional of {Va} at all space points. 
With such relations the macroscopic balance equations (14.121) - (14.141) become hydrodynamic 
equations, i.e. 

%,(r,t) = iV,(r,t I {yj). (5.2) 

The conditions for the existence of constitutive equations then constitute the context for a 
hydrodynamic description. Certainly, it cannot be expected in general that the fluxes can 
be characterized entirely by the hydrodynamic fields for all length and time scales. On the 
other hand important examples exist, such as the pressure tensor for an equilibrium fluid in 
an arbitrary external potential (density functional theory). 

The connection of this question to the statistical mechanical basis for hydrodynamics fol- 
lows from the fact that the cooling rate and fluxes are linear functionals of the solution to the 
Liouville equation, i.e. averages of the form (14. 2p . Therefore, a sufficient condition for con- 
stitutive equations is for the distribution function to be characterized by the hydrodynamic 
fields. A class of "normal" distributions is defined by the functional forms 

p„(r,t) = pjr| {yj) (5.3) 

This means that all time dependence and all the breaking of translational invariance occurs 
only through the hydrodynamic fields. A familiar example of a normal distribution for real 
fluids is the local canonical distribution 



Pe£ (r I {Va}) = exp |g - y drT ^ (r, t) (e (r) - /x (r, t) n(r)) | 



(5.4) 



where g is a normalization constant, and /i(r, t) is the chemical potential (as a specified 
function of the density and temperature for the hydrodynamic fields chosen here). If a normal 
solution to the Liouville can be found then the constitutive equations follow immediately 

C(r, t I {y^}) = 3^^^ ^^y^^ j dTpn (r | {?/„}) w{r) (5.5) 

q(r, t I J) = j dVp^ (r I J) s'(r) (5.6) 

P„^(r, t I J) = j dTp^ (r I {y J) h'^piv) (5.7) 

The origin of hydrodynamics has now been "reduced" to finding conditions for the ex- 
istence of a normal solution to the Liouville equation. Its time derivative in the Liouville 
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equation can be expressed in terms of the hydrodynamic equations (I5.2p 

Substitution of (15.31) into the Liouville equation gives the form for normal solutions 

/ dr ^^,\M r, t I J) + Ip, = 0. (5.8) 
J dya (r, t) 

This intimate connection of constructing a hydrodynamic description and finding a normal 
solution to the Liouville equation is in fact a single self-consistent problem. For specified 
fields, (15.81) is an equation for the F dependence of the normal phase space density as a 
function of the fields. This dependence then allows determination of the normal forms in 
(15. 5p - (15. 7p . Finally, solution of the hydrodynamic equations (15. 2p . with suitable initial 
and boundary conditions, provides the explicit forms for the fields, and completes the nor- 
mal solution. The existence and determination of this solution is the central problem for 
establishing a hydrodynamic description for both normal and granular fluids. 

The concept of a normal solution and its use in the macroscopic balance equations makes 
no special reference to the possible inelasticity of collisions. Nor does the concept refer to 
states near homogeneity or the requirement of representation as local partial differential 
equations. As will be seen below, the familiar Navier-Stokes equations represent a special 
case of this more general idea. For normal fluids, the simple form of the Navier-Stokes 
equations applies for a wide range of structurally simple fluids, with rheology as a counter 
example for more complex fluids. As noted in the Introduction, even structurally simple 
granular fluids can exhibit behavior like complex real fluids for which the Navier-Stokes 
representation fails [23]. However, the generality of the discussion here shows that the failure 
of a Navier-Stokes approximation should not be taken as the absence of a more complex 
hydrodynamic description. 

To clarify the conditions under which a normal solution could be expected, consider 
first a normal fluid with elastic collisions in an initial non-equilibrium state with specified 
hydrodynamic fields {ya (r, t = 0)}, whose values vary smoothly across the system. In each 
small cell the phase space density p(F, t) approaches a local Gibbs distribution characterized 
by the hydrodynamic fields at its central point r, such as is given by (15. 4p . However, this 
is not a solution to the Liouville equation due to the differences in hydrodynamic fields in 
different cells. The solution has additional fluxes driven by these gradients for subsequent 
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exchange of mass, energy, and momentum to equilibrate these fields to their uniform values 
(or to steady values if the system is driven). The first stage, approach to a universal form 
for the velocity distribution, occurs after a few collisions. This establishes the normal form 
of the solution where the hydrodynamic fields and their gradients characterize the state. 
Deviations from the Gibbs density are due to fiuxes of mass, momentum, and energy across 
the cells. These fiuxes are proportional to the differences in values of the fields (i.e., to 
their spatial gradients). The second stage is the evolution of the distribution through the 
changing values of the fields, according to the hydrodynamic equations. 

This two-stage evolution can be expected for granular fiuids as well. The initial velocity 
relaxation will not approach the local Gibbs density, but some other corresponding local 
normal state determined from the inelastic Liouville equation (see below). Subsequently, 
the deviations from this local normal state characterizing spatial inhomogeneities will again 
be via the macroscopic balance equations for the granular fiuid. This is the space and time 
scale for a hydrodynamic description. 



VI. NAVIER-STOKES APPROXIMATION 



The self-consistent solution to (15.81) and determination of Na{r,t \ {Ha}) is a formidable 
problem in general. Specific cases of interest may provide simplifications that allow further 
progress. Consider the example of uniform shear fiow, where the system is driven by Lees - 
Edwards boundary conditions (simple periodic boundary conditions in the local Lagrangian 
frame 2^). At the macroscopic level this state is characterized by a uniform density and 
temperature, and a constant y derivative of Ux - the shear rate. This is an example for 
which all spatial gradients vanish, except the first order derivative of Ux- The latter can be 
small, so that all properties depend nonlinearly on the shear rate, but clearly the problem 
is considerably simplified. 

The class of states to be considered here are those for which all spatial gradients of first 
order can occur, but which are small and all higher order derivatives are negligible. These 
are weakly inhomogeneous states. It is expected under these conditions that the normal 
solution to the Liouville equation can be represented by an expansion to first order in the 
gradients 

p„ (F I {2/4) = p,, (F I {y J) + (F, r | J) ■ Vy„ (r, t) + - (6.1) 
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The dots denote contributions from 'Vya'Vya, VV?/a and higher order gradients. The 
reference state p;^^ (F | {ya}) is the local homogeneous cooling state, analogous to the local 
equilibrium state for a molecular fluid. As a normal state, it is a functional of the exact 
hydrodynamic fields {i/a} and must yield the exact averages for the associated microscopic 
phase functions 

J dV {p^ - Pf^g) aa = 0, aa = {n{r) ,e{v) ,g{r)) . (6.2) 
For spatially constant fields it must reduce to the HCS of (13.211) 

Ph IZ/O/?}) = Pm (r I {1/0/3 + Oy^]) \5y=o, j: = / rfr — \sy=Q, ■■ 

(yyoa J oya [r) 

(6.3) 

where yo^ denotes arbitrary homogeneous values. Thus, the choice for the local reference 
state is not arbitrary and other perturbations of the HCS consistent with (16. 2p are not 
consistent with the second term of (16.11) being of first order in the gradients. Further 
characterization of the local HCS distribution is given in Appendix O The determination of 
Ga (F, r I {ya}) is fixed by the requirement that (16. ip be a solution to the Liouville equation 
(15. 8p to the same order in the gradients. This is discussed in the next section. 

Small gradients means that the relative change in the hydrodynamic fields over the largest 
microscopic length scale io is small: fof^rlnya << 1. There are two characteristic length 
scales, the mean free path and the grain diameter. For a dilute gas the mean free path is 
largest, while for a dense fluid the grain size is largest. In many cases, the condition for 
small gradients can be determined by control over the initial or boundary conditions and it 
turns out to be the usual experimental condition for simple normal fluids. As noted in the 
introduction, some special states for granular fluids entail a balance of the intrinsic internal 
cooling and the boundary or initial conditions such that control over the gradients is lost 
(two examples are noted in the discussion). Thus a careful analysis of each case is required 
before making the assumption of small gradients. In the following it is assumed this has 
been assured. 

The constitutive equations follow from substitution of (16.10 into (I5.5p - (l5.7p 

C(r, t I J) ^ 0(r, t I J) + 3^(^ ^) / dTw{v)Ga (F | J) • Vy^ (r, t) (6.4) 

q(r, t I J) ^ q,(r, t \ {y J) + j rfFs'(r)G, (F | {^4) ■ Vy^ (r, t) (6.5) 
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P„^(r,t I {yj) ^ Peap{r,t \ {yj) + J dTh'^f,{r)G, {T \ {t/J) • Vy, (r,t) (6.6) 

The subscript i on the first terms denote their local HCS averages, i.e. fl5.5l) - fl5.7p evaluated 
with p„ (r I {ya}) Pm (r I {Va})- By construction, these constitutive relations are correct 
to first order in the gradients. However, the fields at different space points of the functional 
also differ by terms of first order in the gradients. Therefore, when an average at some chosen 
space point is calculated a further simplification is possible by retaining only terms of linear 
order in gradients at that point. To illustrate, consider some local property represented by 
a(r, r). Its local HCS average can be evaluated to first order using 



Ph£ 



= PHi{yM,t)}) + (MA{yair,t)}) -r^P^^^^ (6.7) 
where the dots denote terms of higher order in the gradients, and M,^ is defined by 

= / dr' ( r' (6.8) 

Then the average of a(r, r) becomes 

J dTa{T,r)p^,{T \ {y^{r,t)+Sy^}) ^ J dTA{T)p^{T, {y^{r,t)}) 

+ J dTa{T,0)Mf,-Vy^{r,t) (6.9) 

The first term is the local HCS functional evaluated at the "constant" value of {ya{r,t)}, 
in which case it becomes the HCS cooling function of these values, according to (16.31) . 

(r I {ya{r,t)}) = (r, {ya{r,t)}) . The second term is the contribution from 6ya (r',^) 
which is of first order in the gradient ^y^ (r, ^)- 

The first terms on the right sides of (I6.4p - (I6.6I) therefore have two terms at this order, one 
evaluated at the HCS and the second of first order in the gradient. These latter combine 
with the second terms of (I6.4l) -( l6l6l) . These expressions then have their final forms to first 
order in the gradients 

C(r, t I J) U{ya{r, t)}) + Cuiiyair, t)}) V ■ U(r, t) (6.10) 

q(r,t I {yj) ^ -A({i/„(r,t)}) VT(r,t)-/x({y,(r,t)}) Vn(r,t) (6.11) 
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-viiVair, t)}) (^dM^, t)+dM^, t)~V ■ U(r, t)6^p^ (6.12) 

Fluid symmetry (rotational invariance of p^) has been used to recognize q/t(r,t | {ya}) = 0, 
and that contributions to the scalar ( from gradients of the density and temperature must 
vanish. Similarly, contributions to the vector q from velocity gradients must vanish, and no 
second order tensor contributions to Pa/3 can be constructed from density and temperature 
gradients. Equation ( 16.12^ has the familiar form of Newton's viscosity law, while ( 16. lip is 
Fourier's law with an additional contribution from the density gradient. 

These results, together with the macroscopic balance equations f l4.12l) -( HrT^ provide the 
closed set of hydrodynamic equations with constitutive relations calculated determined to 
first order in the gradients. "Determined" means that formally exact expressions are now 
available for calculation all parameters in these equations. For example, ph and are 
properties of the HCS 

Phiiya}) = ^ j d'^Ph (r; {Va}) H'^^ = j drh'aair). (6.13) 

UiVa}) = j dVp^ (F; {yj) = j drw'ir) (6.14) 

Recall that a prime on the phase function denotes the local rest frame, — > = — 
Us(r, t). Equation ( 16.13P now defines the hydrostatic pressure for a granular fluid as a 
function of the local density and temperature. The scaling property of pf^ displayed in 
f l3.2ip confirm that Phi{ya}) oc T and (f^ oc T Similarly, expressions for the transport 
coefficients (^, A, p, k, and r] as phase space averages follow from this analysis. Their 
explicit expressions are deferred to the next section. 

It is appropriate to register at this point the explicit form for the full nonlinear granular 
Navier - Stokes equations [251] 

A^ + nVr-U = (6.15) 

(A + C.)T + A(.+ ^C.+ Q.-)v.u)v.u 

-^{V {daUf3 + dpUa) daUp + V ■ (AVT + pVn)) = 0, (6.16) 
DtUa + {mn)-^da [p - + V • - {mn)-^dpr^ {d^Up + d^U^) = 0, (6.17) 
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These are almost the same as the Navier-Stokes equations for a molecular fluid, except for 
the presence of the cooling rate C/i and two new transport coefficients, C{/ and /i in the 
temperature equation. Perhaps the most significant of these is the cooling rate, which leads 
to new instabilities and new stationary states. To illustrate, consider small homogeneous 
perturbations {6ya} of the HCS solution. The linear equations for the perturbations are 

dtn = = dt6U^, (^dt + Coh + ^^^) + ^^^^^ = (6-^^) 

Using the scaling Coh oc y/Th, and introducing the dimensionless variables 6T*/6T/Th, 6U* = 
SUa/vh{t), these become linear equations with time independent coefficients that are easily 
solved. It is found that there is one decaying mode, one constant mode, and three growing 
modes. A similar result is found for finite spatial gradients, where the same modes are 
unstable at sufficiently long wavelengths. A second interesting effect of the cooling rate 
is the existence of new steady states, that are possible when external work done on the 
system or energy input is balanced by the inherent cooling from collisions. For example, the 
temperature equation for steady, simple shear flow becomes 

^h = ^vdyU.. (6.19) 

It is noted in the discussion that for many of these steady states, it is not possible to control 
the size of the spatial gradients and higher order hydrodynamic effects beyond Navier-Stokes 
order are required. 

To summarize, a hydrodynamic description for normal and granular fluids has been given 
from the exact macroscopic balance equations and the assumption of a normal solution 
to the Liouville equation. For the class of fluid states with small spatial gradients in the 
hydrodynamic flelds, this normal solutions is constructed to leading order in the form (16. ip 
which provides constitutive relations in terms of the flelds and their gradients. In this way, 
the macroscopic balance equations become a closed set of hydrodynamic equations. The 
parameters of this description (e.g., pressure, transport coefficients) are given in terms of 
averages over the solution ( 16. ip . It remains to make that solution more explicit and this is 
the topic of the next section. 
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VII. CONSTRUCTION OF THE NORMAL SOLUTION AND TRANSPORT CO- 
EFFICIENTS 



The objective now is to construct a solution of the form fl6.1l) that is exact up through 
first order in the gradients. This extends recent results for a corresponding solution to 
the Liouville equation resulting from small initial perturbations of the HCS 
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12 . In 



that case the linearized Navier-Stokes equations are obtained for states close to strictly 
homogeneous systems. Here that assumption holds only locally, so the full nonlinear Navier- 
Stokes equations result. The transport coefficients obtained by both methods are the same, 
except for the values of the fields on which they depend. 

As described above, the ultimate use of this solution is to calculate local properties of 
the form 

A{v I = J dTa{T,v)p^{T \ {y„(r, (t)}) . (7.1) 

Therefore, in the following analysis the gradient expansion are referred to the field point r of 
interest, i/a = ya{j^)+5ya- Of course the results will be general and applicable to any choice 
for r. Consider first a general solution of the form 

p (r, t I {y^ it)}) = p,, (r I {y^ (t)}) + A (r, t I {y^ it)}). (7.2) 

The notation makes explicit the fact that there is both explicit time dependence and that 
which occurs through {ya{r,t)} in both p and A, while p^^ depends on t only through 
{ya (r, ^)} by construction. The Liouville equation gives 

d,A+ [ dr N^ir,t\{y^{t)}) + LA = - [ dr-^^N^{v,t \ {vait)}) -Lp,,. 

(7.3) 

It is understood here that the time derivative is taken at constant {y^ {r,t)}. It is expected 
that on some time scale this time dependence goes to zero and (17.31) becomes the same as 
(15. 8p for a normal solution. 

An approximate solution is sought by expanding in the gradients. The right side of (17. 3p 
is evaluated in Appendix [E] with the result 

dtA+ [ dr ^M r,t I {y^ (t)}) + (l -V)LA = -{1-V) (F, {y^ {r,t)})-Vy, (r,t) 

J dyc[r,t) 

(JA) 
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with "fluxes" X^, given by 



(r, {y^ (r,t)}) = {Ct5u, - K^, {{Va (r,t)}) (r, {y^ (r,t)})) 



M,(r,{|/,(r,t)}) = j dr 



r. 



(7.5) 



^Va (r, y Sy=0 

Also £t is the operator of f l3.27p with Th replaced by T (r, t) and is the transpose of the 
matrix 



K 





dn dT ^ 



\ 



\ 



(7.6) 





The matrix K generates the solution for homogeneous perturbations of the hydrodynamic 
equations 



dt5y,{t) = N,{T,t\{y^{T,t)+Sya{t)})-N,{r,t\yo.{r,t)) 
^ -K,^{{y^{r,t)})6y^{t). 



(7.7) 



Further interpretation of its occurrence here is given below. There are two important obser- 
vations that are a consequence of the special choice of Pf^g as the reference state. First, the 
right side of fl7.4p is proportional to the gradients. This admits solutions of the form (16. ip 



A (F, t I {y^ (t)}) ^ (F, t, {y^{r, t)}) ■ Vy„ (r, t) 



(7.8) 



The second observation is the occurrence of the orthogonal projection (1 —V), deflned below, 
which excludes the invariants of the generator for dynamics in (17.40 . This will be shown to 
be essential for the existence of a normal solution. 

Substitution of the form (17.81) into (I7.4p leads to an equation for Gq, (see Appendix[E]for 
details) 

dtG (t) + (1 - P) {ICt + K^) G (t) = - (1 - V) T, (7.9) 

where a compact matrix notation has been introduced. The time derivative is again taken 
at constant {yo, (r,t)} and the dependence on these flelds has been suppressed since they 
are simply parameters of the equation. Integrating this equation with the choice G (0) = 
(initial local HCS) gives the formal solution 



G (t) = - (1 - P) / dt'e~^"^^^'^^y (1 - V) T. 
Jo 



(7.10) 
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Use has been made of the property 

{1-V)[iZt + K'^)V = (7.11) 

which is also proved in the Appendix. 

To interpret the role of (1 — P) in these expressions, defined in fl7.16p below, consider the 
average of some property represented by X (F) 

(X; t) = {X),, + dt'C^ (t')) ■ Vy. (7.12) 
where the first term is the local HCS average and the correlation function (t) is 

C;f(t) = -y"drx(r)(i-p)(e~('^^+^^)*(i-p)T(r)) . (7.13) 

The only explicit time dependence on the right side of fl7.12p occurs in the upper limit of 
the time integral. Suppose that (t) —* for t » tm- Then for t » tm the upper limit 
of the time integral can be extended to oo and the average of X becomes normal - all of its 
time dependence is through the parameters (i",t)} 

(X; t) ^ (X)„ = {X),, + (^^°° dt'C^ (t')) ■ Vy. (7.14) 

Thus sets the time scale for the normal solution, and hence for the onset of hydrodynamics. 
A necessary condition for the existence of a tm is the convergence of the time integral in 
fl7.14p . This means there should be no invariants of the generator for time dependence, 
ICt + in its domain of action. It is shown in Appendix [E] that such invariants exist 

(ICt + K'^)^ = 0, vl/, = ^^M). (7.15) 

The existence of these invariants of (17. 150 provide a microscopic representation of these long 
wavelength hydrodynamic excitations at long wavelength in the spectrum of Ct [3]. The 
effective generator for the dynamics, Ct + , gives a dynamics with these homogeneous 
perturbations subtracted out. However, the projection operator (1 — V) projects out such 
contributions from the definition ( 1E24I) 

{l-V)X = X-m^^ j dVA^X, i j dVA^^p = S^p. (7.16) 

The biorthogonal set are linear combinations of the total particle number, energy, and 
momentum defined in (1E19P . According to (17. 7p K is the generator for homogeneous hydro- 
dynamic perturbations. 
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A technical complication is the occurrence of period time dependence, the Poincare recur- 
rence time. This can be removed by considering the thermodynamic limit of oo, — *^ 
oo at constant N/V . Therefore, the normal solution to the Liouville equation to first order 
in the gradients is written 

Pn (r I {y. m) = (r | {y^ it)}) - {lira j\t'e-i'~^-+^^y (1 - P) ■ Vy.. (7.17) 

It is understood that the limit is taken in the context of an average like fl7.14p with the 
thermodynamic limit followed by the long time limit, all at constant {Uq} = {ya{r,t)}. 

The constitutive equations now follow directly from (16.41) - (16.61) and fluid symmetry to 
get the final forms (I6.10p - ( 16.12^ with the transport coefficients identified as 

CuiiVa}) = J dVW-M,,, + \\m j'^dt'Ciit') (7.18) 

A({yJ) = ^ / dTS- ■M2 + limjyt'CUt') (7.19) 
f^{{ya}) = ^ j drS--M, + limiyt'C^{t') (7.20) 



3V 

viiya}) 



i J dTH-yM,,y + hm dt'C^ (f) (7.21) 
<{ya}) = dVH-^M,,, + hm j'^ dt'C^ (f) (7.22) 



with the correlation functions 



cc {t') = — %- [ driye-(^-"0*' ( e-(^^-+^")*' (1 - V) r) (7.23) 

3nTV J V / 3xx 

(t') = — - f dTS-( e-(^^^+^^)*' (1 - V) r) (7.24) 
3V J \ J 2 

C^^ it') = --L j dTS ■ (^e-(^^^+^'^)*' (1 - P) t)^ 

C" {t') = [ dVH^y L-i^'^T+^^y (1 - V) t) (7.25) 

V J \ J Zxy 

it') = -^j dVH~,, (e-i'-^-+^^y (1 - V) t)^^^ (7.26) 



and 



j dv{l- Pt) «;(r), S = j dv{l- pt) s'(r), Hij = j dv [l - V^) /i^(r) 

(7.27) 
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The adjoint projection operator V"^ is defined by J dVXVY = J dT (P^X) Y. The first 
terms of fl7.18p - (I7.22p come from the expansion of p^^ to first order in the gradients, (16.71) . 
They vanish for non-singular, conservative forces but are non-zero for granular fluids and for 
the elastic hard sphere fluid. Also, they vanish in the low density limit but can be dominant 
at high densities. 

The appearance of Ct in 07. 9p . fl7.17p and the correlation functions for the transport 
coefficients is awkward as it implies using the temperature as a dynamical variable in ad- 
dition to the phase space variables F. This difficulty can be avoided by transforming to 
the dimensionless representation described at the end of Section IIIII The transformation is 
described in Appendix |E] with the results 

(ds + (1 - V*) (iT - A*)) G* = - (1 - V*) [iT - A*) M* {{y^}) , (7.28) 

where an asterisk denotes the function, operator, or matrix expresses in terms of the dimen- 
sionless variables. The matrix A* is 



A* 



\Ch 



(7.29) 



Vo -icy 

Note that now Ct has been replaced by the phase space operator C of (13.271) and the 
temperature no longer appears in this dimensionless form. The variable s is the same as 
that of (I3.35p . with Th T, and represents the average number of collisions per particle. 
The corresponding solution to the dimensionless Liouville equation is 

P: (s, r I {y: (t)}) = pl, (r* I {y: (.)» - (luaj^ ds' (e-(^^*-^*)^' (1 - V*) T*) J ■ V^yt- 

° (7.30) 
The eigenvalues of A are 0, ~\C*h- These are excitations for small homogeneous per- 
turbations of the hydrodynamic equations in this representation. The invariants of (I7.15p 
become 

(^IC* - A*) ^* = 0. (7.31) 

The interpretation of (1 —V*) excluding these invariants of the dynamics is the same as 
discussed above. For practical purposes, both in theoretical and simulation applications, it 
is usually most convenient to use the dimensionless forms. 
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This completes the formal derivation of the nonlinear Navier-Stokes equations, including 
expressions for the pressure tensor, cooling rate, and energy flux including contributions 
up through first order in the gradients of the hydrodynamic fields. These expressions are 
functions of the hydrodynamic fields to be determined by their detailed calculation. The 
scaling property of hard spheres allows the temperature dependence to be obtained directly 
from dimensional analysis, but the density dependence requires confrontation of the full 
many-body problem for calculation. 

A. Example: Shear Viscosity 

To illustrate the results for the transport coefficients, the shear viscosity is considered in 
more detail. The dimensionless form is used, but to simplify the notation the asterisk is 
suppressed. The shear viscosity is found to be 

v{{yo]) = ^J dVH^yM^ '^'"^V I J dTH,ye-(^-^'>'Tr,. (7.32) 
with A3 = — |C/t 

1 ^ 

M, = --J2 qxsd^^sPh, ^v={^- A3) M,. (7.33) 

^ s=l 

The projection operators vanish in this case from fluid symmetry. The volume integrated 
momentum flux H^y is given from flB14p as 



N ^ N 

Hxy = ^ VrxVry ~ 2 ^ QrsxT{r, s)Vsy. (7.34) 



r=l 

In the elastic limit 

N 

2 



1 ^ 

{C - A3) ^ L, p„ ^ -H^lyP^, Mr,^ QxsVys, (7.35) 



s=l 



This result for the flux follows from ( ]A21l) 



T, L (M,pJ = (M,Lp,) + (L_M,) = -H^yP^. (7.36) 

Interestingly, this introduces the generator for time reversed dynamics and the associated 
momentum flux -ff^Tj/- The expression for the shear viscosity becomes 

viiVa}) j dTH^yM^^y + hm 1 ds' j dTH^yC-^^' R-y. (7.37) 
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This is the Green-Kubo shear viscosity for a normal fluid. It differs from that for a fluid 
whose particles interact via non-singular forces by the presence of the flrst term on the right 
side, and by the difference between the two fluxes (forward and reversed dynamics) in the 
flux-flux time correlation function. As noted above, this is a peculiarity of non-conservative 
and/or singular forces. 

Comparison of (17.321) and (I7.37P shows the most signiflcant new affects of dissipative 
dynamics: 1) one of the fluxes, is generated from the reference homogeneous state. This 
representation holds as well for normal fluids, but for granular fluids the homogeneous state 
is no longer given by an equilibrium distribution. Hence the form of the flux is quite different 
from that for a normal fluid; 2) the generator of the dynamics, (£ — A3) , is modifled in two 
important ways. The operator C generates the usual hard sphere trajectories, but modifled 
by an additional scaling operator to account for the dominant effects of coUisional cooling. In 
addition, the dynamics of homogeneous perturbations of the reference state are compensated 
by the subtraction of an appropriate eigenvalue for that dynamics, A3. This effectively shifts 
the spectrum of (£ — A3) to exclude the homogeneous dynamics of the reference state - 
cooling and its homogeneous perturbations. 

B. Impurity Diffusion 

Perhaps the simplest example of hydrodynamics is the diffusion of an impurity particle in 
a host fluid. To illustrate briefly, the host fluid is taken to be in its HCS and unperturbed by 
a single impurity particle. The impurity is taken to be a hard sphere, although its mass, size 
and restitution coefficient for collisions with particles of the host fluid may be different from 
those for fluid particle pairs. The macroscopic balance equation is that for conservation of 
the probability density P(r, t) for the location of the impurity at time t 



This follows from averaging the corresponding microscopic conservation law, where P{r,t) 
and the flux J(r, t) are 




(7.38) 





(7.39) 



32 



The subscript zero distinguishes the position and velocity of the impurity particle. The 
average flux J(r, t) is given by Pick's law in the hydrodynamic limit and for small gradients 



J(r,t)^-D({yJ)VP(r,t). 



(7.40) 



The above analysis can be applied directly to this case as well with the resulting Green- 
Kubo expression for the dimensionless diffusion coefficient 13] (asterisks suppressed) 

mTo 



D = lim 



2moT JO 



ds'C,, 



(7.41) 



The normalized velocity autocorrelation function is 
^ / N _ (vo(s) • vo) _ 2moT 



dV [e' vo) ■ vop;,. 



(7.42) 



This expression for the diffusion coefficient is very similar to the corresponding result for 
normal fluids except for the changes noted above for the shear viscosity - a different generator 
for the dynamics and a different reference fluid state. In addition, the appearance of two 
different temperatures here is due to the different mechanical properties of the host and 
impurity particles (e.g., mass, size, restitution coefficient). This is a reflection of the failure 
of equipartition for granular mixtures, as is expected for any non equilibrium state 26j. 

A simple estimate for the diffusion coefficient can be obtained from truncation of the 
cumulant expansion for the correlation function pj|, I27l | 



exp 



oo 

y- 



Up [s) 



lp=i 



-LUlS 



(7.43) 



where the first cumulant is 



< (Cvq) ■ Vo > 



ly 1 + 



mTn 



1/2 



(7.44) 



{vD 2^" ■ ' V' ' ^oTj 

The arrow indicates an evaluation of the average using the first term in an expansion of the 
reduced pair distribution function for the impurity and host particle in Sonine polynomials, 
and a neglect of velocity correlations. The dimensionless collision frequency u and cooling 
rate evaluated in the same approximation are 

2 21/2, 



27r(l + ao) 



m 



3r(3/2) m + mo 



a 



TX 



3r(3/2) 



X(l-a' 



(7.45) 
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and Xo is the pair correlation function at contact. This gives the diffusion coefficient as 



t^Tq I ^ f ^ mT( 



It is seen that tm = cu^^ sets the time scale for cross over to a normal solution and hydro- 
dynamics. 

This result has been compared with molecular dynamics simulations for the special case of 
self- diffusion (mechanically equivalent particles) 13| . It is found to give excellent agreement 
at low densities and weak inelasticity, while differences grow at both higher densities and 
stronger inelasticity. This result also agrees with the Enskog kinetic theory for both normal 
and granular gases. 



1/2 ^ "1 



VIII. KINETIC THEORY AND BOLTZMANN LIMIT 

The above analysis has provided an exact formal derivation of the hydrodynamic equa- 
tions and associated expressions for the transport coefficients in terms of HCS correlation 
functions. The difficult many-body problem is postponed to the final stage of evaluating 
these expressions. Approximations are introduced only at this final stage. An alternative 
approach is that of kinetic theory where the many-body problem is confronted at the outset, 
and approximations introduced at an early stage. Of course, the two should yield equivalent 
results to the extent that consistent approximations are used in each. In this section, the 
ideas behind a kinetic theory derivation of the hydrodynamic equations are reviewed and an 
approximation expected to be valid for low density gases is introduced. Although attention 
is focused here on low density, kinetic theory methods with different approximations can 



be applied to higher density fluids as well; for references to granular gases see [28|. The 
macroscopic balance equations are obtained from the resulting kinetic equation, and the 
associated constitutive equations are expressed in terms of the solution to that equation. 
The construction of a normal solution is described in analogy to that given above for the 
Liouville equation, and expressions for the transport coefficients of a low density gas are 
obtained. 

The motivation for a kinetic theory representation is based on the fact that hydrodynamic 
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fields are averages of sums of single particle functions of the form 

N 

A = ^a{xr), (8.1) 

r=l 

where Xr < — > (qr, v^) is used to denote both the position and velocity of particle r. The 
particle average of such functions can be reduced to an average over the single particle 
phase space 

{A;t) = J dxia{xi)N j dx2..dxNp(^,t) = n J dxia{xi)f^'^\xi,t), (8.2) 

where f^^\xi,t) is the first member of a family of reduced distribution functions obtained 
by integrating out some degrees of freedom 

( ) N\ f 

ri''f^'"''-'{xi,..,Xm,s)= ^j^_^y_ J dxm+i--dxNpiT,s). (8.3) 

Clearly, f^^\xi,t) is proportional to the exact probability density to find a particle with 
(qi,vi) at time t, regardless of the phase of all other particles. The representation (18. 2p 
expresses the fact the much less information is required for such averages than is contained in 
the full particle state of the system p(r, t), and that it is sufficient to know the dynamics 
in the single particle phase space. 

Equations for the dynamics in the reduced distribution functions follow from the Liouville 
equation by integrating it over some degrees of freedom 



a* + 5^ V, ■ V. - Tih j) ] f^^^Kxu ..,x^, t) 

i=l i<j 

m 

n 

1=1 



m „ 

J dx^+iT(2,m + l)/(™+i)(a;i,..,x™+i,t). (8.4) 



The left side of this equation describes the dynamics of m particles, just as the Liouville 
equation for an isolated system of m particles. The right side represents the effects due to 
interactions with the remaining N — m particles. It describes this interaction for each of 
the m particles with another with Xm+i times the probability density that there is such a 
particle with this phase. The latter is the joint probability for m + 1 particles. It is seen 
that the equation for f^"^^ is coupled in this way to that for This family of equations 



is known as the Born, Bogoliubov, Green, Kirkwood, Yvon (BBGKY) hierarchy 10|]. The 
simplicity of the representation (18.21) is somewhat misleading, since the many-body problem 
is simply transferred to the problem of solving this hierarchy. 
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There have been many attempts to obtain approximate solutions for f^^\xi,t). Most 
of these entail a "closure" approximation so that the first hierarchy equation becomes a 
closed equation for f^^^ alone. This is accomplished by some form of functional assumption 
in which it is assumed that the two particle distribution function can be expressed as a 
functional of the one particle distribution 

f^\xi,xi,t) ^ F{xi,X2,t I f'\t)) (8.5) 

If the functional form F{xi,X2,t \ ■) can be discovered then the first BBGKY hierarchy 
becomes a kinetic equation 

(5t + vi ■ Vi) f^'\x,,t) = J{x,,t I /«(t)), (8.6) 

where the collision operator J(xi,t \ f^^\t)) is identified as 

J(xi,t I /W(t)) =n [ dx2T{l,2)F{x,,X2,t \ /«(t)). (8.7) 



In fact, fl8.5p is so general that such a construction is always possible in principle. A 
meaningful physical approximation requires further guidance. Bogoliubov argued 30| that 
in many cases there is an initial "synchronization" time after which this functional exists as 
a time independent functional 

/(2)(xi,xi,t)^F(xi,X2|/(^)(i))- (8.8) 

In this case the kinetic equation is Markovian. Bogoliubov then went on to construct this 
functional formally as an expansion in the reduced density as a small parameter. This 
procedure was formalized by others via cluster expansions, which allowed a more penetrating 
analysis of terms in the expansion. It was found that the expansion is not uniform in 
time, due to an unexpected class of recollisions of particles ("rings") leading to secular 
contributions. A similar detailed analysis for granular systems has not yet been carried out. 

Here, a formal solution to the hierarchy will be constructed for low density gases using 
again the reduced density as a small parameter and for the special case of inelastic hard 
spheres. This expansion is initiated by returning to the dimensionless form of the Liouville 
equation (13.401) . The corresponding dimensionless form of the hierarchy is found to be 



(a. + £;(e)) r^"'\xi,..,xi,s) = Y, / rfwir(z,m + i)r(™+i)(xt 

i=i ^ 



,..,x:^„s). (8.9) 
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with the definitions 



i=l ^ ^ i<i 



(8.10) 
(8.11) 



- n^,J] 



j da e(g*,. ■ ^)(g^ ■ a) [a-2(5(q;. - eo-)^-.^ - 5(q^ + 



cr 



(8.12) 



(8.13) 



The length scale has been chosen to be the mean free path ^ = l/(ncr^), where n is the 
density. Low density is now defined by e << 1, i.e. the grain size is small compared to the 
mean free path. This suggests looking for a solution to the hierarchy as an expansion in this 
small parameter 

= + ^ (8.14) 

It is readily shown 



31| that the hierarchy is solved exactly to zeroth order in e by 



r(m) / 

Jo l-^l) ■ ■ ■ ; ^my ^ 



(8.15) 



1=1 



where fQ^^\xi, s) is the solution to 



ids + vt . V,*) f;^'\xl, s) + ^V.* ■ [v*/*(^)(a;t, .)J = J*K I /o^'^ (s))- 

The operator J* is the Boltzmann collision operator for inelastic hard spheres [iC 

I /o^'^ is)) = I dx;T;{l,2)C\xl,s)C\x;,s), 

where Tq(1,2) is T*{i,j) at e = 

T;{z,j) = j da e(g., ■ ^)(g., • 6-)5(q*^.) \cr%l - l] . 

The effect of different locations of the colliding particles (contained in the Boltzmann- 
Bogoliubov collision operator SOj) does not contribute at this order. Finally, C,*Qh is given by 
fl3.3ip at e = 0, which reduces to (using the symmetry of the HCS) 



(8.16) 



(8.17) 



(8.18) 



a = (1 - ^ / dvldvlf;'^\vl)f;^^\vl)gll 



(8.19) 



37 



where f^j^^ is the stationary HCS solution to the Boltzmann equation (18.161) 



2 



V*/, 



i J Oh 



(8.20) 



This provides a complete description for the determination of fQ^^\xl, s) and the calculation 
of averages such as (18.21) at low density. It is an exceptional non-trivial result, and somewhat 
special to the hard sphere interactions. A similar analysis for non-singular forces leads to a 
collisionless mean field theory. 

It is sufficient for the purposes here to stop at this point and inquire about solutions 
to this kinetic equation leading to the hydrodynamics appropriate for a low density gas. 
However, it is instructive to digress briefly for a discussion of higher order terms. These 
consist of two types, those that correct the single particle distribution in the product form 
(I8.15p . and those that generate correlations. One class of terms of the first type are those 
coming from the e dependence of T*{i, j) which give the "collision transfer" corrections to the 
Boltzmann equation. This suggests a "renormalization" of the expansion to one in powers of 

at constant T*{i,j). Then, to lowest order the above results are obtained again but with 
the replacement Tq(1,2) T*{i,j). This is the Boltzmann-Bogoliubov low density kinetic 
theory that retains the different locations of the colliding particles. The first corrections to 
this theory have the form 



/r(™)(xt, .., s) = j2i[ f<'\x:, s)f<'\x;, s) + J2i[ f;^'\xi s)g%xi x*, .), (8.21) 

j=l i^j i<j kj^ij 

where the expression for f*^"^^ holds for m > 2. Thus, the reduced distribution functions for 
any number of particles are determined as a sum of products of the single particle functions 
fo^^\^i^ s) and fi^^\xl, s), and pair function (j*(x*, 0:2, s). These are determined from the 
set of three fundamental kinetic equations [31] 



(a. + vt.V,.)/o*(^'(xt,.) + %V.* 



Xi, s) 



dx;r{i,2)ff\xi,s)f, 



X 



21 • 



{ds + CI) fl 



(1), 



dxlT\l,2)G\x\,xl,s). 



{d, + CI + CI) G*{xl, xl s) = r (1, 2)/o«(a;t, .)/o* 



f*(i), 



(8.22) 

(8.23) 
(8.24) 
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The operator is defined over functions of xl by fl8.22l) 

dx;r (1,2) {^C\xi,s)h*ix;) + h*{xi)f;^'\x;,s)) . (8.25) 

These low density results f l8.22p -( l8.24p are remarkably rich. The first equation for the 
single particle distribution function is given by the Boltzmann equation, as expected. The 
second equation provides corrections to the Boltzmann results in terms of the two particle 
correlations. This provides a means to understand the limitations of the Boltzmann equation. 
The third equation gives the dynamics of these pair correlation. They are driven by a source 
term which is the result of the binary collision operator acting on uncorrelated single particle 
distributions that are solutions to the Boltzmann equation. In more detail, it is possible to 
show that the solution to this pair correlation equation, together with the equation for f^^^^ 
give corrections to the Boltzmann equation due to correlated recollisions ("ring" collisions) 
among many particles. These effects dominate at long times and therefore show that the 
Boltzmann theory is not accurate at asymptotically long times, no matter how low the 
density. The third equation also gives the means to study the dynamics of pair correlations, 
which also are characterized by correlation recollisions. Finally, not shown is the equation 



for pair correlations at two different times, which also follows from this analysis 31 1. 



A. Hydrodynamics from Kinetic Theory 



Consider now the Boltzmann kinetic equation for the one particle distribution function 
f l8.22p . The dimensionless hydrodynamic fields, scaled relative to their values in the corre- 
sponding HCS, ^M), are 



n*(r*,s)T*(r*,s) 
n*('r*,s)U*('r*,sl 



I 1 \ 



1*2 
3^ 



r(r*,v*,s) 



(8.26) 



Multiplying by fl8.22p by 1, f *^,and v* and integrating, successively, gives the macroscopic 
balance equations. They are the same as those of fl4.12p - fl4.14p . in dimensionless form, 
except for expressions defining the cooling rate, heat flux, and pressure tensor. Here they 
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are integrals over /*(a;*, s) 

C* = (1 - I dvldvinv\ V*, s)r{v\ V*, (8.27) 

P*. = 2 j dv*V*V* f*{r*,Y*,s), q* = j dY*V*^V*f*{r*,^r*,s), (8.28) 

and V* = V* - U*(r*,s). 

The conceptual discussion in the sections above showing that hydrodynamics follows from 
the existence of a normal state applies here as well [29] . A normal solution to the Boltzmann 
equation for first order in the gradients (Navier-Stokes hydrodynamics) is obtained in a way 
analogous to that for the Liouville equation in Section IVIII [32!] • The solution is written as 

r (r*, V*, .) = f;,{V*, K(r*, .)}) + A(r*, v* | {yHs)}). (8.29) 

The reference state //^ is the local HCS distribution given by 

fkiV*, K(r*, .)}) = n*ir*, s)f*{ J^ ), (8.30) 

A/i *(r*, s) 

where fl{v*) is the HCS solution given by (18.201) . Because of this choice for the reference 
state it is found that A(r*, v* | {?/q(s)}) is of first order in the gradients 



{K + 62uCT* + w* -Vr^yl) (8.31) 



dy* 

The subscript on Cqi means that the spatial gradient contribution to in (I8.25P is 
excluded {C^-^ = Cl — v* ■ Vr*)- the first equality, use has been made of the fact //^ is 
related to the HCS solution by (18.301) . The second term contains the hydrodynamic fluxes 
N* + 62uCT* (recall equation (15. 2p ) which are first order in the gradients. As in Section IViTt 
it is essential that the reference state is an exact solution to zeroth order in the gradients to 
assure that A is of first order 

A(r*, V* I K(.)}) ^ G:(., V*, K(r*, .)}) ■ V*y: {r*, s) . (8.32) 

Substituting this form into (18.301) . the analysis proceeds in a similar manner to that described 
for the Liouville equation in Appendix [El The details will not be given and only the result 
is quoted for the normal solution 

r(v^{y:}) = /4(l^^K})- (^iim^'rfs'(e-(^^o-^*>'(i-pr)T*)J ■v*y:. (8.33) 
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The corresponding normal average for some property X (v) 

{X; {y:}): = {X; {y:})l + (^^" dsCf {s; K})) ■ Vyt, (8.34) 

Cf {s; {y:}) = -ld^rX (v) (1 - VI) (e-(^^o-^*)^ (1 - VI) T* (v))^ . (8.35) 

These results should be compared with those of Section IVlIl to note that they have the same 
form as those for the Liouville equation. The single particle projection operator now projects 
orthogonal to the invariants of [ICqi — A* ), where A* is the same matrix as in (I7.29p . 



0, 



(v) = v]/, (V, {y J) J dY.Ap (Vi, {yj) X (vi) , 

/ n \ 



iT, r\r t l^ 9fih{V,{ya}) 7 r\r r \\ 

(^ya 



3n* 



{V*^ - |T*) 



Also, the flux T* (v) is generated from the HCS in much the same way 



;i - p*) T* (v) 



(8.36) 



(8.37) 



(8.38) 



-(l-n*)(^/A-A*)qvl/„ 
_(l_p*)V^„ 



(8.39) 



The constitutive equations for Navier-Stokes hydrodynamics also are the same as those 
of Section [VIIl and follow from insertion of (18.331) into (18.271) and (I8.28p . The transport 
coefficients are of course different, due to the limitations imposed by the low density ap- 
proximations. For example, the shear viscosity is now 



V*{{ya}) = -lim J ds I dvH^yC' 



-A3 



T 



(8.40) 

r = v^d^J*, (V*) (8.41) 
This is quite similar to the general result (17.321) . The first term on the right side of that 



result vanishes in the low density limit and does not appear here. It has been shown that 
the general result reduces to this obtained from kinetic theory when the correlation func- 
tion there is evaluated with the same low density assunaptions as applied here 33|]. The 



details for other transport coefficients are described in 3^, l35|] and will not be repeated 



here. It is found that the results agree in detail with those obtained from the alternative 
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Chapman- Enskog method to construct a normal solution (properly adapted to the granular 



Boltzmann equation 36|]). The Chapman- Enskog representation requires solution to an in- 
tegral equation. This has been done approximately, using a truncated polynomial expansion 
for the solution (or equivalently, using the cumulant expansion approximation described for 
the diffusion coefficient above). The Green-Kubo representation given here provides another 



approach to evaluation using direct Monte Carlo simulation of the correlation function 



IX. DISCUSSION 



3- 



The presentation given here has addressed the origins of hydrodynamics for a granular 
fluid. There have been three components to the discussion: 1) an exact set of balance equa- 
tions, based on the starting microscopic description - either Liouville dynamics or kinetic 
theory; 2) the concept of a normal state which implies constitutive equations in terms of 
the hydrodynamic fields, and hence a hydrodynamic description; and 3) an approximate 
construction of the normal solution for small local gradients, leading to the Navier-Stokes 
hydrodynamic equations. Only the last component entails specific limitations on the state 
conditions considered, with explicit neglect of higher order gradients in the fields. It is 
presumed that these state conditions can be controlled to assure small gradients, and then 
Navier-Stokes hydrodynamics would apply. Below, it will be noted that some new kinds 
of steady states for granular fluids preclude the control over gradients and that conditions 
for Navier-Stokes hydrodynamics cannot be attained in these cases. But there are many 
other cases, verified in both simulations and experiments, for which a Navier-Stokes order 
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4l| . It is worth empha- 



hydrodynamics is found to give an accurate description [38, 
sizing again that the failure of a Navier-Stokes description does not mean the absence of a 
hydrodynamic description, only that the constitutive equations are more complex. 

Even when the gradients are small, it is expected that the hydrodynamic description 
will dominate all other "microscopic" excitations only on some sufficiently long time scale. 
Thus the characteristic hydrodynamic times should be long compared to the decay times for 
all other excitations. For normal fluids this is assured by the fact that the hydrodynamic 
fields are associated with global conserved quantities (number, energy, and momentum), 
whose lifetimes become infinite in the long wavelength limit, while all other characteristic 
times remain finite. This is not the case for granular fluids, since the energy is no longer 
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conserved. Instead, in the long wavelength limit there is a new characteristic time deter- 
mined by homogeneous cooling, (^^^. Some have argued that for this reason the energy or 
temperature should not be considered among the hydrodynamic fields for a granular gas. 
However, there are many good reasons why the temperature should remain one of the fields 
in the long time, long wavelength description of granular fiuids. For example, consider the 
simplest case of a general homogeneous state. On the time scale of a few collisions it ap- 
proaches the HCS for which the temperature obeys Haff's law (13.341) . The latter is the exact 
hydrodynamic description in this case which is seen to dominate even though it has a time 
scale Ch^- For spatially inhomogeneous states, it has been shown here that the dominant 
background cooling can be suppressed by the dimensionless representation. This means that 
the relevant approach to hydrodynamics is relative to this background state. The picture 
described above of each cell in the fiuid rapidly approaching a local HCS due to velocity 
relaxation, followed by approach to spatial uniformity by fiuxes between the homogeneously 
cooling cells is consistent with this. 

A more precise study of this problem is possible within the simplifying features of the 
Boltzmann kinetic theory. There the response to an initial small spatial perturbation of 
the HCS can be studied to see if there is a dominant set of collective modes at long times 
and long wavelengths. For a normal fiuid this is established by showing that the spectrum 
of the generator for dynamics in the linearized Boltzmann equation has five isolated points 
that are smaller in magnitude than all other spectral points, and therefore dominate at long 
times. These are shown to exist and to be the same as the eigenvalues of the linearized 
hydrodynamic equations. A similar analysis of the linearized granular Boltzmann operator 
has shown the existence of hydrodynamic eigenvalues in its spectrum in the long wavelength 



limit 
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42]. Further progress can be made using kinetic models for this Boltzmann oper- 



ator, showing that the hydrodynamic modes of the spectrum remain isolated and smallest, 
even for strong dissipation 42|, |43|| . In summary, there is very good evidence for the dom- 
inance of a hydrodynamic description for granular gases in the long time, long wavelength 
limit. 

As with normal fiuids, there are states for which the simple Navier-Stokes form for hydro- 
dynamics does not apply because the gradients are no longer small. In the case of normal 
fluids steady states can generally be studied in the Navier-Stokes domain by external control 
of boundary conditions. For granular fluids, there are new types of steady states associated 
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with an autonomous balance of external constraints and the internal cooling. For example, 
a fluid under shear has viscous heating due to the work done on it at the boundaries. Nor- 
mally this is compensated by a temperature gradient that induces a compensating heat flux 
to produce the steady state. However, a granular fluid can compensate for the external work 
done via its collisional cooling. The latter scales as T^/^ so any amount of work done can 
be accommodated by the system choosing an appropriate steady state temperature. This 
balance is given by the energy balance equation (14.131) whose steady state form becomes, for 
uniform temperature and density 

C(n, Ts)Ts = -^PapdaUp. (9.1) 

This imposes a relationship of the given boundary shear, the coefficient of restitution a, and 
the steady state temperature T^. Any attempt to produce small gradients by decreasing the 
boundary shear also decreases Ts at constant a. Since the dimensionless measure of small 
shear scales as T~^/^ it is found that this dimensionless shear can never be brought within 
the accuracy of the Navier-Stokes hydrodynamics [44]. Another example of this is a granular 
fluid pinned between two walls at constant temperature. The internal cooling introduces a 
temperature gradient toward lower temperatures between the walls. The gradients estab- 
lished are controlled by the cooling and cannot be made small by changing the boundary 
temperatures. It can be shown that the steady state temperature proflle is never that of the 



4J]. Such discrepancies have been seen in recent 



Navier-Stokes hydrodynamic prediction 
molecular dynamics simulations 45j. 

These last observations justify characterizing granular fluids as "complex". Under some 
conditions they are well-described by the local partial differential equations of Navier-Stokes 
form; under others they have a rheology or complex dissipation that is not observed in 
normal fluids, or only for those with structural complexity. Studies to date suggest that a 
hydrodynamic description in its most general sense - closed equations for the hydrodynamic 
flelds - is a reasonable expectation, although the associated constitutive equations may be 
complicated and state dependent. Discovery of generic constitutive equations (e.g., for shear 
flow) is both the challenge and opportunity posed by granular fluids. 
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APPENDIX A: HARD SPHERE DYNAMICS 
1. Generator of Trajectories 

The phase function A (F) in (13. 4p is evaluated at a phase point = 
{qi(t), . . . , qAr(t), vi(t), . . . , vjv(t)} which has evolved from its initial point at t = 0, 
F = {qi, . . . , qTv, Vi, . . . , vjv}. Therefore A (Fj) can be considered as a function of the 
initial point and time, A (Fj) = A {T,t). The simplest example to illustrate this is a gas of 
non interacting particles for which Ft = {qi + Vit, . . . , qAr + v^t, Vi, . . . , vtv}. This corre- 
sponds to free streaming of all particles at constant velocities. A phase function at time t 
can therefore be expressed in terms of a generator for the dynamics in the following way 

A (Ft) = A ({qi + vit, ...,qN + VArt, Vi, . . . , vn}) 

= e^^'A ({qi, . . . , q^, vi, . . . , v^}) , (Al) 

where L is the sum of generators for translations for each particle coordinate 

^o = J]v,-Vq,. (A2) 

r 

In this form the dependence on F and t is made explicit. 

In the presence of interactions among the particles the time dependence of the positions 
and velocities is more complex, but still given by a deterministic rule for evolution of the 
initial point. If the particle interaction is pairwise additive and due to conservative, non- 
singular forces F(gjj), then (lAip is replaced by 

A(Ft)=e^*A(F), (A3) 
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where the generator L 

1 ^ 

L = Lo + - $^m-^F(g,,) ■ (V^. - V^,) . (A4) 

It is easily verified that this generator changes the positions and velocities according to 
Hamilton's equations {L is the linear operator representing the Poisson bracket of its 
operand with the Hamiltonian for the system). Generalization to include non-singular, 
non-conservative forces is straightforward. 

Hard sphere dynamics is also represented by a rule for evolution of the initial phase point: 
particles freely stream until one pair is at contact. At that time the velocities of that pair 
are changed according to a specified rule (e.g., that of (13.21) ) constrained to conserve total 
momentum for the pair. Subsequently, all particles freely stream again until another pair is 
at contact, when the corresponding velocity transformation of that pair is introduced. The 
generator for this process cannot be represented simply as in ( ]A4[) using the corresponding 
singular force. To discover the correct form [21!], assume the first collision is between particles 
1 and 2, and let Ti(qi2,gi2) denote the time at which it occurs as a function of the initial 
separation and relative velocity. Therefore, the time evolution of any phase function A (F) 
can be given compactly by 

A (Fi) = (1 - (t - n)) A (Fi) + e (t - n) A (F,) + - (A5) 

where 9 is the Heaviside step function, 

/ 1, X > 1 \ 

e(x)= - . 

\ 0, X <l 

The dots in (lASP denote contributions that arise only on the time scale of the third and later 
collisions. The first term on the right side of the first line contributes only for times before 
the first coUision, so that qi,2(^) = qi,2 + ^i,2t, ^i,2{t) = ^i,2- The second term contributes 
only after that collision, so qi,2(i) = qi,2 + ^1,2^1 + 2 (^ ~ '^1) > ^i,2(t) = v'^ 2- The prime 
indicates that the velocities of the colliding pair have been changed for t > ti according to 
the collision rule (13. 2p 

v'l = vi - i (1 + a) (qi2 • gi2) qi2, V2 = V2 + ^ (1 + a) (qi2 ■ gi2) qi2- (A6) 

This distinction between F = {qi(t), . . . , qAr(t), vi(t), . . . , VAr(t)} for t < Ti and t > Ti 
refiects the discontinuity in the velocity at t = ri. 
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Differentiation of (1A6P with respect to time gives 

dtA (r,) = (vi(t) ■ Vq,(t) + V2(t) • Vq,(i)) A (1*) + 5 {t - T,) (AA (F.j - A (r.,_,)) + - 

(A7) 

The first term is defined for times before and after the colhsion, with appropriate changes in 
qi^2(^) and Vi 2(t) as given above. The second term occurs only at the time of the collision 
and is proportional to the discontinuity in A A (Ft-J = lime_»o {A (Ft-J — A {Tt^^^)) due to 
the instantaneous change in the velocities of the colliding pair. This can be represented by a 
substitution operator 612 (gi2(Ti)) as in (13. 9p . but here defined for functions of the position 
and velocity at time ri 

(F,J = (612 (gi2(ri)) -I) A (F.,_,) . (A8) 

Note that 612 (gi2(T)) changes Vi^2(^) at constant qi^2(^)- The time dependence of the delta 
function can be expressed through qi2(i) — qi2(''"i) = gi2(i — ^i), where gi2 is the relative 
velocity before the collision. By definition of the hard sphere collision qi2(Ti) = cr is a vector 
of length a implying contact between the two particles. This can occur only if the particles 
are directed toward each other, i.e. for a'-g,i2 < 0- Thus the delta function in time can be 
written in terms of a delta function for position 

V |cr-gi2| / 

= (-^•gi2) 1^ ■ gi2| 5 (gi2(t) - a) . (A9) 



Hence ( 1A11I) takes the form 



d,A{r,) = (vi(t)-Vq,(,)+V2(t)-Vq,(i))A(F,) 

+S (quit) - ff) (-S--gi2 (n)) |ct • gi2 (ri)| 

X (612 (gi2(ri)) -1)A{T (n)) + - (AlO) 

Since the delta function enforces t = ti the latter can be expressed equivalently as t in the 
second term. With the representation flASP the generator of trajectories L is identified as 



Vl-Vq,+V2-Vq,+T(l,2) + -- 

1 ^ 

Lo + 7;J2'^ir,s) (All) 



2 

r,! 
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with the binary coUision operator T{r, s) given by 

T(r, S) = 6 {qrs - Cr) {-^rs ■ Zrs) |^rs " g,rs\ (Ks - 1) • (A12) 

The second hne of (I A lip extends the analysis to longer times where there is sufficient time 
for many collisions. 

2. Two adjoint generators 

Next, consider the generator L defined in (13.101) . This can be identified from 
J dTA{T)Lp{T) = - j dV {LA{T)) pT) 

= dvp (r) (v, ■ Vq,A (r)) r^P (r) s i^rs - a) 

r=l r,s 
xQ{-ars ■ grs) \^rs " grs | [Ks " 1) ^ (r) . (Al3) 

Define the inverse of brs by h^}hrs = KsKs = 

KsSrs = g'rs = Srs - (1 + «) ft"^ {^rs " grs) ^rs- (Al4) 

A useful identity is given by 

j dTX (T) hrsY (r) = j dVa-^X (fe^/r) Y (F) . (A15) 

This follows by changing integration variables from (v.^, v^) to (fersV,,, hrs^s)- The factor 
is the Jacobian for this change of variables. Also 

{^rs ■ grs) = -a'^ars " grs- (A16) 

With these results f lXTSl) becomes 



j dTA{T)Lp{v) = fl[J rfrp(r)(v,- Vq,A(r)) + ^rfq,- J rfv,rfr,^,(v,p(r) A(r)) 

+ I dTA (r) 6 {qrs - (t) l^rs ■ grs I 

r,s 

X (e {ars ■ grs) - e {-^rs " grs)) P (F) . (A17) 

An integration by parts has been performed and the formal adjoint Liouville operator is 
identified as 

1 ^ - 

L = Lo--Y,T{r,s), (A18) 



2 

r,s 
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where, the new binary colhsion operator is 

T(r, S) = 6{qrs - Cr)\grs ■ q.rs\ (© (qrs " grs) Ci'\} " © (-Qrs " grs)) • (A19) 



A second adjoint generator can be interpreted as that for reversed dynamics [12] • It is 
defined by flSTHjl 

e-~^ (p(r)B(r)) = (e-*V(r)) (e-*^-5(r)) , (A20) 

or equivalently 

I {piv)Biv)) = (Ip(r)) BiT) + p(r) (L_i?(r)) , (A21) 



The distributive property of flA22p is clearly satisfied by Lq, so L_ can be written in the 
form 

1 ^ 

L_^Lo--Y,T4r,s). (A22) 

■r,s=l 

The new binary operator T-{i,j) is identified from from the corresponding condition 

piT)T4r,s)BiT) = T(r, (p(r)B(r)) - (T(r, .)p(r)) s(r) 

= S{qr,s - cr)a"^67/|g,s ■ qrs|0(-grs • qrs) (p(r)5(r)) 
-fi(r)(5(g,rs - o-)a"^67/|g,.s ■ qrs|6(-grs ■ q,.s)p(r). 

= [SiQrs - a)a''b;}\Srs ■ q.s|e(-g,, ■ q..)p(r)] {b;} - l) B{T) 

= piT)6iqrs - a)|g,, ■ q..|e(g,, • q,,) {b;,' - l) B{T). (A23) 

The last equality follows from a boundary condition for hard particles at contact 
S{qrs - (r)a~'^b~^\grs ■ qrs|0(-g„ ■ qrs)p(r) = S{qrs - cr)\grs ■ qrs|0(grs ■ qrs)p(r). (A24) 

This can be understood from the fact that the fiux of particles moving towards each other 
at contact must equal the fiux of particles separating at contact. Since the velocities of the 



former and latter are, respectively, g^^ = &j/gij and gij, this gives 46|] 

S{qrs - (T)\g',, ■ qr.|e(-g;, ■ q,,)p(r')c/r' = S{qrs - (T)\grs " qr«|e(g,, ■ qrs)p{T)dT. (A25) 
This result implies ( ]A24l) . The form for T-{i,j) is then identified from ( 1A23I) as 

T_(r, S) = 5{qrs - (y)\grs ■ q,rs\Q{grs " ^.rs) {b~^ - l) (A26) 
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APPENDIX B: BALANCE EQUATIONS AND FLUXES 

In this appendix, the forms of the fluxes and sources associated with the local number 
density n, energy density e, and momentum density g are identified from the microscopic 
balance equations associated with these densities. Consider first the case for phase functions 
of the form 

N 

A{r, t) = e^*A(r), A{r) = J] a(v,)<5(r - q,). (Bl) 

s=l 

Its forward time evolution is 

{dt-L)A{r,t) = (B2) 
The action of the Liouville operator is 

N 



LA{r) = -Vr- Xlv.o(v.)5(r-q,) 

1 ^ 

+ - J2 s) (5(r - q,)a(v,) + S{r - q,)a(v,)) (B3) 

r,s 

Next, use the identity 

d 

S{r-c[r) = 5(r-q5)+ / d'j—S{r-cis + ^cisr) 

Jo ^1 

= 5(r-q,) + Vr^ • / d7(5(r-q, + 7q,r)q,r, (B4) 
Jo 

to obtain the microscopic balance equation for A{y, t) 

5tA(r,t) + Vr-B(r,t) =5(r). (B5) 
The flux B(r, is identified as 

B(r) = Vfca(vjfc)5(r - q^) - - / d'j^S {r - qr + TQrs) <lrsT{r, s)a{vs). (B6) 

fe=l "^0 r,s 

The first term on the right side is the "kinetic" contribution due to the kinetic motion of the 
particles, while the second term is the "coUisional transfer" part of the flux that requires no 
motion. The source term in the balance equation is 

TV N 

2 



TV N 

-^(r) = 9 E E - "Irmr, s) (aK) + a(v,)) (B7) 



'■=1 j¥=s 
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It vanishes for any two particle "collisional invariant" a(vs) +a(v^.) which is the same before 
and after a binary coUision. For elastic collisions or conservative forces these are given by 
a(v,) = 1, v,t;2. 

Consider now the local microscopic number, energy, and momentum densities defined by 

N N ^ N 

n = ^(5(r-q,), e = ^ -mv^5(r - q,), g = ^ mv,5(r - q,). (B8) 

s=l s=l s=l 

The microscopic balance equations become 

dtn{r, t) + m-^Vr ■ g(r,t) = (B9) 

dte{r, t) + Vr ■ s(r,t) = w{T,t) (BIO) 

dtgair,t) + Vr,Kp{v,t) = (Bll) 

The absence of sources on the right sides of ( 1B9I) and ( IBllI) occurs since 1 and v are still 
summational invariants for inelastic collisions. The loss function in the energy equation is 

N 



1 /I 1 \ 

r,s ^ ^ 

1 ^ 1 

= -J2S{r-ci^T{r,s)-mgl. (B12) 

r,s 

The second equality follows from the fact that the center of mass energy is conserved, by 
momentum conservation. The energy and momentum fluxes are 

^i^) = ^^s-mvl6{r - qs) - - ^7 ^ 5 (r - q^ + 7q^,) q„T(r, s)-mv^. (B13) 

s=l r,s 

hapir) = ^ mVsaVsf^Sir ~ ^r) - ^ j ^7 ^ 5 (r - q^ + 7qrs) qrsaT{r, s)mVsf3. (B14) 

s=l ''^ r,s 

1. Time reversed dynamics 

Now consider the same phase function as in ( IBll) except for the reversed trajectory (still 
with t > 0) 

(-9t-L_)A(r,-t) = 0, v4(r,-t) = e-^-M(r). (B15) 



N N 



Ev..V,,-i5^T_(r,s) (B16) 



2 

s=l r. 
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The derivation of balance equations for the backward dynamics is similar to that above with 
the result 

- dtA{v, -t) + Vr ■ B_(r, -t) = S.{v). (B17) 

The flux B_(r) is 

N , „1 TV 



N ^ „1 1\ 

B(r) = ^ v,a(vs)5(r ~ ^s) - ^ ^7 ^ 5 (r - + 'yqrs) <irsT-{r, s)a(v^ 

.8 = 1 "^0 



(B18) 



and the source 5'_(r) is 

N N 

2 



1 r,s 

Interestingly, the fluxes and sources for the forward and backward balance equations are not 
the same due to the differences between T(r, s) and T_(r, s). This is a peculiarity of hard 
spheres for both normal and granular fluids. It occurs for non-singular, non- conservative 
dynamics as well. 

The corresponding microscopic balance equations for reversed dynamics are 

9n(r, —t) 



dt 

9e(r, —if:) 
dt 

dg^iv, -t) 



+ m-^Vr • g(r, -t) = (B20) 

+ Vr-s_(r,-t) = ?i;_(r,-t) (B21) 
+ Vr,/i-a/3(r,-t) = (B22) 



dt 

1 ^ /I 1 \ 

^-(^) = "2 5Z -^(^ - ^) + ^rnvlj . (B23) 



The energy and momentum fluxes are 
s_(r' 



1 1 1 

) = X] v,-mv^5(r - ^s) + - ^7 ^ 5 (r - q,. + 7qr,) q,.,T_(r, s)-mv2. (B24) 

TV 1 /-l ^ 

h^ai3{r) = ^ mt;^at;s/35(r - q^) + - / ^7 ^ 5 (r - q,. + 7q.rs) qrsaT-{r, s)mVsp. (B25) 

Note that the fluxes are not the same as those for the forward time conservation laws, (1B13P 
and (lB14p . because of the differences between the binary collision operators T(r, s) and 



T_(r, s). This observation is not purely of academic interest, as both sets of fluxes occur 
in the Green-Kubo expressions for transport coefficients for a normal fluid (see for example 
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APPENDIX C: LOCAL REST FRAME TRANSFORMATION 



The microscopic densities and fluxes can be represented in terms of contributions from 
local convective flow plus their values in the local rest frame. This is accomplished by a 
partial Galilean transformation on all velocities, = — U(r, t), where U(r, t) is the 
macroscopic flow field defined by Eq. (14 .Op and the point r is the same as that of the local 
property being considered. A straightforward calculation for the densities gives the results 
n(r) = n'(r), Z(r) = /'(r) and 

e(r) = e'(r) + g'(r) • U(r, t) + ^mn{r)U\r, t), g(r) = g'(r) + mn(r)U(r, t). (CI) 

A prime on a phase function indicates that same function evaluated at ^ v^. The energy 
and momentum fluxes transform according to 

s(r) = s'(r) + U(r,t) (^e'(r) + g'(r)-U(r,t)+^mr2'(r)f/2(r,t)^ 

+ \g'{r)U'{r,t) + h'^,{r) (C2) 

V(r) = h'^pir) + g'M)Up{r, t) + %ir)Uo.{r, t) + mn'(r)f/„(r, t)f/^(r, t) (C3) 
The averages of these are then n{r,t) = n'{r,t), {w{r);t) = {w'{r);t) and 

e{r,t) = e'{r,t)+^mn{r,t)U\r,t), g(r, t) = mn(r, t)U(r, t). (C4) 
(s(r);t) = (s'(r);t)+U(r,t) (^e'{v,t)+^mnir,t)U'ir,t)^ + {h',piry,t) (C5) 

APPENDIX D: LOCAL HOMOGENEOUS COOLING STATE 

The HCS solution to the Liouville equation corresponds to spatially constant hydrody- 
namic fields. The local HCS is a reference state distribution that is not a solution to the 
Liouville equation but approximates the true nonequilibrium normal solution. Physically, 
it represents a partitioning of the system into cells such that each cell is in a HCS at its 
own values for the fields. These values are chosen to be same as the exact values for the 
nonequilibrium state. This condition is expressed in (16. 2p . Recall that the HCS distribution 
has the scaling property expressed by (I3.2ip 

p, {T;t) = (Iv, {t))-'"^pl (1^, ^^^^} , v^{t) = V2T,(t)/m (Dl) 
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This suggests that the local HCS distribution should be defined as 

v{qr,t) = y/2T{qr,t)/m. (D3) 

The dimensionless function pjj has its arguments changed to refiect the fact that the hy- 
drodynamic fields at the positions of each particle can be different. Thus translational 
invariance is broken only by the spatial dependence of the hydrodynamic fields. It is verified 
that fl3.2ip is satisfied, and clearly p^^ — > p^^ if the value of the fields is the same everywhere. 
The stronger requirements in (16.31) are also satisfied. To illustrate, let the reference value for 
the fields be evaluated at a chosen point r and write i/a (qr, t) = ya (r, t) +Sya (qr, t)- Then 
direct calculation gives 

In the same way higher functional derivatives have similar properties 



Vr /JJ'' \sy=o= ^TT^^- (D4) 



6Ui{r',t)6U,{r",t) "^=° dUi {r,t) dU, {r,t)- ^ ^ 

These requirements exclude other choices for p^^ that are simply perturbations of Pf^. The 
essential role of flD4p is demonstrated in the next appendix. 

The modified scaling of Pf^^ (T, \ {ya}) provides alternative forms for the functional deriva- 
tives with respect to temperature and velocity fields. For example, 

SPhe (r \{yf3}) I _ ^ 5\nv{qs,t) 

— ^fj^^ I 5y=o - -2^^^Tp-^Vv, ■ (V, - U (q„t))p;,(|?/^(r,t)|) 

-— J2 ^i^s - r') Vv. ■ (v, - U (r', t)) p, {{y^ (r, t)}) (D6) 



2T(r',t) 



s=l 



Similarly, the corresponding velocity derivative can be written 

SUdr',t) - ^ ^(^^ - ' ) d^, 

These functional derivatives are therefore phase space densities derived from derivatives of 
Pfj^ with respect to the density. 

As stated above, the local HCS distribution is not a solution to the Liouville equation. 
Instead, it differs from a solution by terms proportional to gradients in the hydrodynamic 
fields. 
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APPENDIX E: SOLUTION TO THE LIOUVILLE EQUATION 



In this Appendix the normal solution to the Liouville equation p„ is obtained to first 
order in the gradients. First write it as a deviation from the local HCS defined in Appendix 

P = Pm + A- (El) 
Substitution into the Liouville equation gives (17.31) 

dt/\+ [ dr iV,(r,t| fa^(t)})+lA = - / dr-^^iV„(r,t| {|/„(t)})-Ip,,. 

(E2) 

where the time derivative is taken at constant {y/?}. The assumption is that p^^^ has been 
chosen as the correct reference state in order that A is of first order in the gradients. This 
requires that the right side of (15. 8p should be proportional to the gradients. To show this, 
represent p^^ as an expansion to first order in the gradients about the reference values 
{ya{r,t)} 

Phi = Phi{yair,t)}) + j dr'(^^-^^^ (^^(r',t)_y^(r,t)) + .. 

= pj{y„(r,t)}) + m^(r,{y,(r,t)})- Vy;3(r,t) + - (E3) 

where the dots denote terms of higher order in the gradients, and m,^ is defined by 

^Pm 



m,(r,{y,(r,t)}) = j dr" 



fr"-r) 



With this expansion for p^^^ the functional derivative on the right side can be evaluated 

/ *'r%7T'^"(-'.* I fa"}) = ^^^#T^iV„(r,i I {,„}) 

J 6ya (r', t) dy^ (r, t) 

gmff (r, {ya(r,t)}) . , , r „ , 
+ iV„(r,t I {,.})■ V,,(r,t) 

+m„(r,{y,(r,t)})- ViV,(r,t | {y^}) (E5) 



The last two terms are of first order in the gradients so it is sufficient to use the lowest order 
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form A^„(r, t \ {ya}) -Sa2Ch {{Vair, t)}) T (r, t) to get 



-C,(K(r,t)})T(r,t) 



dva.fi {Y,{y air, t)}) 
dT{r,t) 



The right side of ( lE2p is now 



(r, *)}) ■ "'f '"f " ■ Vy, (r. <) (E6) 

<9l//3 (r, t) 



t^r'^^^^iV,(r',t I {yj) - Lp^, = -/:rp, - {Crnifs (r, {y„(r,t)}) 

+1112 (r, {2/a(r, i)}) a.. ^\ ■ ^^/3 1^' ^) 



^Ph({z/a(r,t)}) 



(iV„(r,t I {y,}) 



dVa (r, t) 
+<5.2C/.({l/a(r,t)})T(r,t)). 



(E7) 



The operator Ct is the same as that introduced in Section IIIIl except with the HCS values 
{Vha} replaced by the true values at the point of interest, {ya{r,t)} 



CtX {{y^ir, t)}) = -C, {{y.ir, t)}) T (r, t) + lX. 



(E8) 



dT (r, t) 

The first term of (]E7p vanishes, CTPh = 0. This is the first important consequence of the 
choice Pf^^ for the reference state; it is a solution to the Liouville equation to first order in 
the gradients. 

Next consider the last term of flE7p and recall that Na{r,t \ {ya}) arises from averaging 
the microscopic conservation laws 



Nair,t I {ya})= dT{Laa{r,{yair,t)}))p^ 



dVaa (r, {ya{r, t)}) iCrpn + Ch ({?/«}) T 



dVaa (r, {y„(r,t)}) Lp^ 
dpn 



(E9) 



where (r, {ya{r,t)}) are linear combinations of the local densities of number, energy, and 
momentum 

n{r) — n (r, t) 
3^(e'(r)-|r(r,t)n(r)) • (ElO) 



/ 



fla (r, {ya{r,t)}) 



\ 



1 g'(r) 



n(r,t)m 
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Note that 



dTac,{r,{ya{r,t)})Ch{{ya})T^ = Chi{y»})T-^ J c/Fa^ (r, {?/„(r, t)}) p„ 

da^ (r, {ya{r,t)}) 

Pn 



<Hi{ya})T J dT 



dT 

6a2Chi{ya})T, (Ell) 



SO 



Na{r,t I {y^}) + 6^2Chi{ya})T = - j dTa^{r,{y^{r,t)})CT{Phi + ^) 
= -j drac,(r, {?/o(r,t)})/:Tm/3(r, {?/o(r,t)}) ■ V?//3(r,t) 

- J dra,(r,{|/„(r,t)})£rA (E12) 



With these results the Liouville equation flE2p becomes 

dtA+ [ dr ^ ,M r,t\ {ya}) + LA - f dra„ (r, {2/„(r, t)}) Z^A 

J Sya{r,t) dya{r,t) J 

^rm^ r, {y^ r, t } + m2 r, {y^ r, t } — -— ■ Vyp r, t) 

oyp (r, t) J 

+^^^^^^iy^ / "^^^^ ^^^^ '^^"'^ ^^^^ ■ ^^'^ ^^^^^ 

The right side is now exphcitly proportional to the gradients, so a solution of the form (16. ip 
is possible 

A = y rfrG„ (r, r, t \ {y^{v, t)}) ■ Vy„ (r, t) . (E14) 

This is a direct consequence of the choice of p^^ as the reference state, as seen by the above 
analysis where all terms of zeroth order in the gradients cancel. Substitution of flE14l) into 
(]E13p . evaluating the functional derivatives on the left side, and equating coefficients of the 
gradients gives the desired equation for (r,r, t | {?/«}) (here and below it is understood 
that all fields are evaluated at y^ = yai^^, t)) 

dtA + (ctG^ (t, {yj) + G2 (t, {ya}) ^^^^^ ^ 

dph {{ya}) 



dTaa (r, {ya}) CtG/s (t, 



dya 

Crnif^ (r, {?/„}) + ma (r, 



dyf3 

dPh {{ya}) 



dyo 



j dTaa{r,{ya}) CTmf3{r,{ya}) (E15) 
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The structure of this equation can be exposed further by introducing a set of functions 
{ipo} that are biorthogonal to the set {da} 

The biorthogonahty property is 



(E16) 



dTa^ {r, {ya})il> air', {y a} [ 



S (r' - r) 6af3 



5yp (r' t) J ^''^ i{y<^})^ 

5dair,{ya}) ,r 

d^— — ^ . Ph [{ya}) 



(E17) 



Integrating over r' gives the related orthogonahty condition 

1 



V 



(E18) 



where 



dya 



{{ya}) = / dr 



N -N 



\ 



1 O/ 



nm 



\ 



J 



(E19) 



Several useful identities identities follow from the condition that the averages of fields 
Ha (r, {Ua}) vanish for p^g^ and therefore, according to fl6.2p 

dTaa{v,{ya{r,t)})p^,{{ya{v,t)}) = Q = j dri, ({y„(r, t)}) p^, ({|/,(r, t)}) . (E20) 



Direct calculation gives the same result for {{Va}) 



dTaa{v,{ya{r,t)})pf^{{ya}) = ^= / dTAa{{ya{r,t)})pf^{{ya{v,t)]) 



This implies 



(E21) 



dVaa (r, {ya{r, t)}) rap (r, {?/^(r, t)]) = ^ = j dTAa ({l/a(r, t)}) (r, {ya{Y, t)}) 

(E22) 

dTaa{v,{ya{r,t)})Gp{t,{ya{v,t)}) = {) = j dTAa{{ya{r,t)})Gp{t,{ya{r,t)}) (E23) 
Finally, define the projection operator V by 

VX = ^Ha^ I dVAaX. (E24) 
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Then (lElSP takes the final form 

{dt + il-V) {ICt + i^^)) G (t, {y^}) = -{l-V) {ICt + K^) M ({yj) (E25) 
Here / is the unit matrix and is the transpose of the matrix 



K 











\ 



dn 









dT 





(E26) 



The matrix K generates the solution to the hydro dynamic equations in the absence of 
gradients, Eqs. (16.181) . 

The explicit dependence on r has cancelled in flE25P (as it must, for a normal solution) 
as a consequence of 

{ICt + K'^) = r {IZt + K'^)^ = 0. (E27) 
The second equality follows from the fact that \1/ form the null space for {ICt + K^) 



[ICt + K'^)^ = 0. 



(E28) 



This can be demonstrated by direct calculation 



dya 



9^2 T 



(E29) 



Equation ( ]E25p is the primary result of this Appendix, giving the exact solution to the Li- 
ouville equation up to first order in the gradients. As a final simplification, a transformation 
from the operator Ct to the phase space operator C of (I3.27P can be made by introducing 
dimensionless variables 

Gi(r,t,{yj) = ^(MT))-'^Gt(r,.,n£=^), 

-3N 



G2(r,t,{yJ) = -{£v{T)r'^G;{T\s,nf) 



G3(r,t,{yj) 



{iv{T)y^^ G; {T*,s,nf) 



rV2 



(E30) 



The first factor in each case gives the dimensions to compensate for the associate gra- 
dient multiplying Gq,. The second time arises in each case because the solution to the 
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Liouville equation is a density in phase space. The dimensionless phase point T* = 
{qt, q^, VJ, .., V^} is defined by 

q: = qr/i, V; = (v, - U(r, t)) (E31) 

where Vh{t) is the thermal velocity defined in (13.221) . with Th T(r, t). With these defini- 
tions 

TV 

v^itr ' 2 



{dt Ir,{,.} -UTdr \r,t) = ^ ( -^,d, |r.,{,.} -ClTdr \r^,t +f E + " ^v,) 



r=l 



i \Vh{t) ' " OT 

N 



-ClTdT \r^,s +y E (3 + " Vv?) j (E32) 

Here s = s{t,T) is a dimensionless ttime scale, and C/l is the constant dimensionless cooling 
rate of fl3.3ip . A judicious choice for the dimensionless time is seen to be 

^ ds 

-dt\r,iy^}-aTj^\td.\r'^d,\r* (E33) 



or 



- ^rft. (E34) 



It can be shown that this definition agrees with that of (13.351) in the sense s{t) = s{t,T{t)) 
dsit) = '-I^dt, m = T(0) (l + ^^t) " . (E35) 



Equation (]E32p simplifies to 



N 



{dt |r,{,.} -UTOt \r,t) = ^ Ir- -ClTdr Ir-,. +y E + " ^v,) j • (E36) 

Now, taking into account the forms (lE30p . the equation (]E25|) for Gq, has the corresponding 
dimensionless form 

(ds + (1 - V*) [iT - A* + K^*)) G* = - (1 - V*) [iT - A* + K^*) M* {{y^}) , 

(E37) 

where an asterisk denotes the function, operator, or matrix in terms of the dimensionless 
variables. The matrix A* arises from the first factors of (]E30P associated with the dimensions 
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of the respective gradients. 







A* = c;; 

ICy 



(E38) 



The notation is simphfied by introducing the matrix A* 







\ 



A* ^K^*-A* = ICh 







(E39) 
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